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ABSTRACT 


In majiy applications involving wave propagation, piolDlem domains are often 
very large or unbounded. A common numerical method used to solve such problems 
is to truncate the domain via artificial boundaries to form a hnite computational 
domain. To accomplish this, Non-Reflecting Boundary Conditions (NRBC’s) which 
minimize spurious wave reflections are imposed. The quality of the solution strongly 
depends on the properties of both the NE^C and the wave behavior. 

This dissertation explores the use of Higdon NRBC’s to solve shallow water 
equations (SWE’s) in a dispersive environment. A linearized SW’E model is developed 
that includes stratification and advectdon effects. Initially a single NRBC is used to 
truncate a semi-infinite channel. Later four NRBC’s are used to restrict an infinite 
plane. In both cases finite rectangular domains are formed, A scheme de\’eloped by 
Neta and Givoh is used to rapidly discretize high-order Higdon NRBC’s. Finite differ¬ 
ence methods and are used in all numerical schemes, which are solved explicitly when 
possible. Results will show that Higdon NRBC’s can be used effectively to restrict 
large rectangular domains when solving SWE’s that include the before mentioned 
effects. 
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I. 


INTRODUCTION 


Phenomena involving the pi operation of waves in imhounded (oi veiy lajge) 
domains are applicable to many fields including acoustics, electromagnetics, meteo- 
lology, and geophysics, Howevei, it is infeasible to compute nnmeiical solutions foi 
legions of this scope, Theiefoie, it is necessary to define artificial boundaries that 
reduce the size of the domain. To accurately model the wave action in the truncated 
region, artificial boundary conditions must be imposed that allow waves propagating 
inside the region to pass freely without spurious reflections, which would pollute the 
computational domain. Such a boundary condition is known as the Non-Reflecting 
Boundary Condition (NRBC) and is the main subject of this dissertation. 

In general, it is not possible to construct a boundary condition that will accom¬ 
plish the criteria perfectly, but during the last 25 years research has been conducted 
to develop NRBC’s that after discretization lead to stable, accurate, efficient and 
easily-implemented schemes [Ref, 1, 2, 3, 4], Investigations in the late 70’s to early 
80’s produced a number of low-order local NRBC’s, e,g. the Engquist-Majda [Ref, 
5] and Bayliss-Turkel[Ref. 6] boundary conditions. The exact non-local Dirichlet-to- 
Neumann (DtN) NRBC [Ref. 7, 8] and the Perfectly Matched Layer [Ref. 9] boundary 
conditions were developed in the late 80’s and early 90’s. Subsequently, higher order 
NRBC’s were introduced, but were difficult to employ beyond the 2^^ or 3^'^ order. 

Only since the mid 90’s have practical higher order schemes been developed. 
Collin 0 [Ref. 10] proposed such a scheme for two-dimensional time-dependent wave 
in a rectangular domain. Grote and Keller [Ref. 11] extended the domain to three 
dimensions in a scheme based on spherical harmonic transformations. They extended 
their work to include elastic waves [Ref. 12]. These findings were independently pub¬ 
lished by Sofronov [Ref. 13] in Russian literature. Hagstrom and Hariharan [Rjef. 14] 
constructed high-order NRBC’s for two- and three-dimensional domains based on the 
analytic series representation for the outgoing solutions of these equations. Guddati 


1 



and Tassonlas [FLef. 15] devised a high-oidei NRBC scheme foi time-depen dent waves 
in a 2-dimensional wave guide using rational approximation and recursive continued 
fractions. Givoh [Rjef. 16] derived high-order NRBC’sfor ageneral class of wave prob¬ 
lems leading to a symmetric finite element formulation. These early investigations 
utilized either time-harmonic waves or non-dispersive time-depen dent waves. 

Wave dispersion, however, is an ever present phenomenon. In the late 80’s and 
early 90’s, Higdon developed NRBC’s for non-dispersive waves [Ref. 17, 18, 19, 20], 
but later showed that his schemes could be applied to the dispersive (Klein-Gordon) 
wave equation [Ref. 21]. Higdon’s work involves low order formulation of his scheme. 
Givoh and Neta [Ref. 22, 23, 24] present an algorithm for implementing the Higdon 
NRBC to any order using high-order FD discretization. They further developed 
methods to rewrite the Higdon NRBC without using high order derivatives and to 
generate Higdon parameters that maximize the non-reflection property of the NRBC 
in a dispersive wave environment. 

In the present work I will develop high order Higdon NRBC schemes for use 
with hnearized shallow water equations (SWE’s) in Cartesian coordinates. The SWE 
model is further enhanced to include the effects of stratification and advection. A 
single Higdon NRBC is initially apphed as an artihcial bcnindary on one side of 
a semi-infinite channel. Later four Higdon NRBC’s are apphed to the sides of a 
rectangular domain to restrict an infinite plane. Finite-difference schemes are used to 
numerically solve the problems. Discrete forms of the Higdon NRBC, based on the 
work of Givoh and Neta, are then employed on the artificial boundary. The results of 
several numerical examples are reported to vahdate the SWE models as weh as the 
use of the Higdon NRBC as an effective means of restricting a very large domain. 
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II. 


MODELING GEOPHYSICAL FLUID 
FLOW 

In this diaprtei we develop a model that desciihes large scale geophysical flow. 
Many of the details of this derivation are developed hy Pedlosky [Rjef. 25]. We 
start hy considering the dynamics of a shallow rotating fluid layer. The fundamental 
condition that characterizes a shallow layer is: 



where D and L characterize the scale of vertical and horizontal motions. This charac¬ 
terization is applicable to large scale atmospheric and oceanic flows where the vertical 
scale of the fluid layer is of the order of miles, while the horizontal scale is of the order 
of hundreds or thousands of miles. Work hy Rjosshy [Rjef. 30] and Stommel [Rjef. 31] 
show that such a working geophysical model can he developed hy assuming that the 
fluid is; 

• Incompressible (density independent of pressure), 

• Inviscid (no internal frictional forces), 

• Homogeneous (not stratihed with regards to density). 

Later in this paper, the “homogenous^ assumption is relaxed. The following physical 
laws are applied in deriving a geophysical model: 

• Conservation of Mass, 

• Conservation of Momentum (Newton’s 2^^ Law), 

• Conservation of Energy, 

• Second Law of Thermodynamics. 

Since the fluid is incompressible, the energy equations are uncoupled from the model. 
Hence, we consider only the conservation of mass and momentum. The basic form 
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□f tliGS0 laws apply to iixGd quantitiGS of mattGi. In thG analysis of Unid flow, wg 
aiG concGinGd with hxGd volnniGS (G.g. thG propGitiGS of a fluid within a control 
volnniG rathGT than propGitiGS of individual particlGS in motion). What follows is 
a dGiivation of thG consGivation of mass and momGntum Gqnations as thGy apply to 
control volnmGS. 

A. CONSERVATION OF MASS EQUATIONS FOR 
FLUIDS IN A CONTROL VOLUME 

ThG consGrvation of mass law statGS that thG total changG of mass in a control 
volnmG must he equal to the net flow of mass Gntering and leaving at the volume 
surface. If /? is the densily of the fluid and pu ■ n is the mass flux at a point on the 
volume surface, then by involdng the conservation of mass law we can write; 

j pdV = - j [pu ■ n)dA, (II.l) 

where V is the volume and ^ is the volume surface area. By the divergence theorem 
we know that: 

/ {fm- n)dA = / V ■ {pu)dV. (IL2) 

Applying this to (II.l) yields: 

j [^ + S7-{pS)\dV = a, (IL3) 

which imphes: 

^ + V-(/«:) = 0. (IL4) 

Since the fluid is assumed to be incompressible and homogeneous (e.g. p is constant), 
this equation is rewritten as; 

V-f? = 0, (II.5) 

or, in three dimensional Cartesian coordinates; 

(II.S) 


du dv dw 
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wheie T£, 1 ', ajid w are the velocity components in the z-, y-, and ^-diiections lespec- 
tively. This expiession is called the contimiity equation. It desciihes the conseivation 
of mass foi a volume element containing an incompressible fluid and is the £ist con¬ 
trolling equation of the geophysical fluid flow model. 

B. THE MOMENTUM EQUATIONS FOR FLUIDS IN A 
CONTROL VOLUME 

The next set of controlling equations is the momentum equation and is based 
on Newton’s Second Law: 


Force = (mas s)(a€celeratio7i). 


Acceleration components must be carefully derived, because geophysical phenomenom 
occur in a non-inertial rotating frame. In the following analysis, we first consider a 
derivation of these equations for a fluid in an inertial two-dimensional Cartesian sys¬ 
tem. A rotational component is then introduced to generate non-inertial momentum 
terms. A third dimension is added and acceleration components in vector form are 
obtained, Finally an Earth model is developed on which the geophysical fluid flow 
equations of momentum will be based. 


1, Acceleration Components for a Fluid in a Cartesian 
Inertial Frame 

We first consider acceleration in an inertial frame (eg. one that is not accel¬ 
erating or rotating). The acceleration component in the x-direction is: 

dv. 

a. = -, ai.7) 


where v. is the velocity in the x-direction. Since we do not treat fluids as individual 
particles, but rather as a continuum of matter, u not only depends on time, but also 
on the spatial components z, y, and .o. Therefore: 


du du du du 
du = —dt -^dz 


(II.8) 
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and (II.T) becomes: 


dv. du dudx du dy du d^ 

^ dt dt dz dt dy dt dz dt^ 


(II.9) 


wbich is Iewritten as; 


du du du 


(II.10) 


wbeie 1 ' and w are components of velocity in the y- and .r-direction respectively. We 


now define the following special symbol; 


D d d d d 


(nil) 


which is often called the deriz'ntiz'e following a fluid [Ref. 27]. Using this symbol we 


write acceleration as: 


Similarlv; 


or in vector form: 


Dr Dw 

— ~~~ 3*lici Q-j — ~~~ . 

^ Dt Dt 


( 11 . 12 ) 


(11.13) 


(11.14) 


We now consider the effects of the rotation on the acceleration vector. 


2 . Acceleration Components in a Two-Dimension 
Cartesian Rotating Frame 

Rotation is an important factor in dynamics if the time to complete a 
single rotation is on the order of or less than the time talcen by an object/hnid 
field to cover a distance L at a speed U ^e.g. — < 1 where = — J . For example, 
consider a large disc that completes one rotation every 5,000 seconds (A^'). On this 
disc a particle travels 10 kilometers (L) at a speed of 1 meters-second“^ (D), it will 
complete its journey 10^ seconds (ti). Thus — = .5 and we conclude that rotation 
will influence the momentum equations. 
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Figiiie 1. Cajt 0 sdaji Two-Dimeiisdonal Rotating Fianie 

In deiiving the lotational acceleiation teinis foi the siniplified two-diniensdonal 
model, tracldng variables is a challenging task. Figure 1, which superimposes a two- 
dimensional Cartesian rotational frame on an inertial frame at a common origin, 
introduces these variables. Additionally, Table I tabulates the variables with regards 
to their reference frames. 

At time t, the rotating x-axis makes an angfe with the inertial X-axis. 


Table I. Variables Used to Derive Rotational Components of Momentum Equation 



Inertial Frame 

Rotational Frame 

a 


angular rate of rotation 

Unit Vectors 

r, j 

h j 

Coordinates 

X, Y 

y 

Velocity Vector 

U 

u 

Velocity Components 

[/, V 

V., r 

Acceleration Vector 

A 

a 

Acceleration Components 

A, B 

a, b 
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Applying equations of lotation, it follows that; 


i = I cos -|-J sin Qt, 
j = —I sin fit-|-J cos fit, 


(IL15) 


and: 


z = X cos fit -1- V"sin 
y = — _Y sin -1-y cos 
The velocityr vectoi in the lotational frame is; 

dz^ dy^ 


(11.16) 


(1117) 


wheie the hist time deiivative yields the velocity components of the lotational frame; 


dz dX dV . 

V. = — = -cos£7t-l- -sin — QJf sin -I-£7T cos 

dt dt dt 


r = 


dt 


±X . dY 

-sin -1--cos fit — r2_Y ^os — fiT sin fit. 

dt dt 


Similarly, the velodly yectoi in the inertial frame can be expressed as: 

- dX^ dV^ 

^ = it'+IT-’= 

Algebraic and trigonometric manipulation of (11.15) yield: 


(11.13) 


(11.19) 


I = icos — j sin Qt, 

J = isin £7t-1-j cos fit. 

Using these and (11.19) yields: 

U = [ cos -|- ^^sin ] id- [ — sin -1- cosf^i ] j. 
\ dt dt i \ dt dt f 


( 11 . 20 ) 


( 11 . 21 ) 


This expression along with (11.15) and (11.13) reveal the following relationships be¬ 
tween the inertial and rotational velocities: 


U = V. — rty, 

V = 1' -|- Qx. 


( 11 . 22 ) 
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We now consddei the TotataonaJ components of acceleration; 


d'^z- Syr T .r 

Taking the second derivative with respect to time of (11.16) yields: 


a = 


SX ^ 

cos Sit -1- ^-;r sm S7t 


b = 


dS 


d\X 


dt‘ 


dX . ^ dV 
—— sm S7t -I- —— cos Sit 
dt dt 


T 2Q 

— Q? cos fit -|- fiy sin fit), 


dt^ 


sin S7t -1- 


d^Y 

dS 


±X 


dY 


cos S7t — 2S7 —— cos S7t -1- —— sin S7t 


dt 


dt 


— Sl^ sin Sit -1- cos fit). 
We express the inertial acceleration vector as; 


d^X . ^ ^ \ r 

sm fit H—^ 7 ^ cos fit j. 


(11,23) 


(11.24) 


(11.25) 


(11,26) 


This and (11.20) result in; 

- fd^X ^ d^Y . ^\r 

A = —— cos fit -1- sin fit 1 -1- . , ^ _ , 

\dt^ di^ j \ dt^ dt^ 

Based on the preceeding, we generate a final simplifying expression that relates the 
rotational and inertial components of acceleration: 

a = A -\- 2QV — fl^ T, 
b = B-\- 2QU - 9}y, 

or comhining this with (11.22) we have; 


(11,27) 


A = a — 2fli' — fl^j, 
B = b 2flii — fl^^, 


(11.23) 


This equation shows the effects of rotation on the inertial acceleration vector. The first 
contrihution proportional to fl and velodly is called Coriolis acceleration. The second 
contrihution proportional to fl^ and the coordinates is called centrifugal acceleration. 
We now refine these equations for the Earth model. 
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3, Acceleration Components in the Earth Model 

FiQm hencefoith we shall use the teims geophysical and “large-scale^ inter¬ 
changeably, but to be more precise a quantitative definition is offered. Let: 

_27r_ 

tinie to complete one revolution/rotation' 

For example, = 7.27radians-second"^ (one rotation every 24 hours). As 

before, we define; 


L = the length of a phenomenom, 

U = the speed of a particle within the phenomenom. 


Using £7, L and U we introduce a dimensionless parameter called the Rjossby number: 


U 

2aZ’ 


(11.29) 


If the Rjossby number e is of order one or less, the phenomenom is “large-scale^ and 
the Earth’s rotation is a significant factor in the momentum equations [Rjef. 25]. For 
example, a 500 kilometer long ocean current with a speed of 15 meter-second"^ has 
a Rossby number € = .205. Therefore the current is large-scale (or geophysical) and 
is significantly affected by Earth’s rotation. 

In order to model the Earth, the momentum equations for a rotating sphere 
are now considered. This non-inertdal frame is compficated by the fact that at any 
given point on the sphere we perceive ourselves on a planar disc. To recreate this 
“human experience^ we set up a three dimensional Cartesian system whose origin is 
at the point of the observer. The z-, y, and .:r-axis are positively oriented to the east, 
north, and “straight-up^ respectively. The Earth’s axis of rotation however, is none 
of these, but rather passes through the North and South Poles. This conundrum, 
depicted in Figure 2, is the reference frame in which we develop the geophysical fluid 
flow equations. 

Before continuing further, we revisit the two-dimensional rotating system pre¬ 
sented in the previous section in Figure 1. We add a third dimension to this system 
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Figure 2, The Earth Model 


by dehniiig a unit vector k at the origin that is perpendicular to the plane. A rotation 
vector is described as Sl = STk. We use this to write the inertial velocity equations 
(11.22) in vector form: 


u 






= 




V 


1' T S7x 




(11.30) 


Now let r = and consider: 


n X r = 


i j k 

0 0 
x ^ 0 


Q(-^iT xj). 


(11.31) 


This allows ns to write the inertial velocity equations (11.22) in condensed vector 
form: 


U = u T n X r. 


(11.32) 


Similarly, the inertial acceleration equations (II.2S) in vector form are: 



-2S^i' 


X 

A = a -|- 

2S^u 

- 



(11.33) 
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Using an analogous appioacih we show that the condensed vector form of the inertial 
acceleration equations is: 


A = a-1- 2n X u -1- n X (n X r). 


(IL34) 


In this form, the Coriolis acceleration is 2n x u and the centrifugal acceleration is 
n X (n X r). With respect to (11.32) and (11.34), the time derivative for the inertial 
frame is equivalent to applying the operator: 


Dt 


(11.35) 


B 


where — is the derivative following a fluid introduced in (II.11). 

We now consider two additional simplifications in the development of the ac¬ 
celeration equations for the Earth model. First we neglect extraneous terms resulting 
from the Earth’s curvature. In general this can he done if L r where r is the ra¬ 
dius of the sphere. On Earth, L < 1000 kilometers is acceptable [Rjef. 2S]. A second 
simplification to (11.34) comes from our intuition about planetary phenomenom. The 
centrifugal force is an outwardly normal force, or from the vantage point of the Earth 
observer, a force that acts straight up. Yet nothing ever is “hung^ upward from the 
face of the Earth because gravitational forces keep centrifugal forces in check. In the 
absence of rotation, gravity would hold the Earth together as a perfect sphere. The 
presence of rotation and accompanying centrifugal forces distort the sphere, flattening 
it to the extent that gravity and the centrifugal force negate each other. Hence, we 
neglect centrifugal forces in the Earth model and the inertial acceleration equations 
become: 

A = a — 2flr, 

B = b -\- 2r2i^, 

or in vector form: 

A = a -|- 2n X u. (11.37) 

Using these equations, we continue with the Earth model. 


(11.36) 
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Ref 0 TTiiig to FigiiTG 2 we wiite the Totataon vector H in terms of the unit 
vectors i, j, and k: 

n = cos -|-f^sin 0k, (11.38) 

where 0 is the degree of latitude of the observer. Since i is always orthogonal to Sl, 
it does not appear in the equation. Using (11.35) we write; 

Du 


a =--1- 2n X u, 

Dt 


where; 


2n X u = 


i j k 

0 2r2cos0 2Qsin0 

V. ^ w 


or; 


Th erf ore; 


D^ 

0 ^ = - -I- 2S7^i? cos 0 — 2£7r sin 0, 

Dt 


(11.39) 


(11.40) 


2n X u = (2^iw cos 0 — 2^ti- sin 0)i -|- 2^lv. sin 0j — 2^tu cos 0k. (II'^1) 


D^ 

Oy = H- 2Q.TI sin 0, 


(11.42) 


Dw 

= —— — 2i7iicos0. 

Dt ^ 

For convenience, we define the following quantities; 

Coriofis Parameter; / = 2f^sin 0, 

Rjeciprocal Coriolis Parameter; = 2f^cos0, 
allowing us to rewrite (11.42) as; 

_L f f 


(11.43) 


d-tf — 


Dr ^ 


(11.44) 


Dw 

^3 = - ft^ 

^ Dt 


13 



The angle and theiefoie the Coiiohs paranietei, is positive in the Noithein Hemi¬ 
sphere and negative in the Southein Hemispheie. The Coriolis parameter is zero on 
the Equator where ^ = 0. The reciprocal Coriolis parameter is positive everywhere 
except at the poles = ±90) where it is zero. 

One final simplification is now applied. As stated earlier, our model of geo¬ 
physical fluid dynamics is based on the shallow water assumption where the vertical 
dimension (D) is very small relative to the horizontal dimensions (L)\ 

f- 

In other words, geophysical fluid flow is “almost two-dimensional^ and therefore ver¬ 
tical flow components are negligible (e.g. w and si 0). Thus we rewrite 

(11.44) as; 

(11.45) 

These are the acceleration components of the Earth model that are used to generate 
the momentum equations for the geophysical fluid flow model. 

4, Forces Acting on a Fluid Control Volume 

Two types of forces act on a fluid in a control volume: Body forces which are 
proportional to the volume mass and surface forces which are proportional to the 
volume surface area. 

Gravity, the only applicable force in the Earth model, acts along the z-axis 
toward the center of the Earth. It acts on a rectangular volume element with sides 
dz, dy, and (iz and is given by; 

Bgravi-ty = —mg = —{p dz dy dz)g = —pgdV, (11.46) 


Dv. . Di' 

= + a, = 0. 


where dV is the incremental volume of the fluid element. 

Salient surface forces include pressure and viscosity. Pressure acts in a direc¬ 
tion normal to the volume’s surface as shown in Figure 3. The total difference in 
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Figure 3. Piessuie Effects on a Rectajigular ETluid Volume 


piessnre along dz is; 



To obtain the total force due to pressure in the j-direction, this expression is multi¬ 
plied b^ the surface area (dy dz) upon which it acts. The result is; 


pr<=££tir<=o. 


= -^4. 4y d. = 

ax ax 


(11.47) 


This force acts in the direction of the negative gradient (e.g. from an area of high 
pressure to an area of low pressure). Similarly the expressions; 


"pressure^. — 






(11.43) 


describe the force due to pressure in the y- and s-directions; 

Viscosity manifests itself in surfaces forces called shear. Figure 4 summarises 
the viscous stress tensor. Viscous stresses are symmetrical (e.g. = r^^). A summa¬ 

tion of viscous forces in the j-direction vields; 


xy 






dy T 


dr^ 


d- 


■d^ dx dy, 


(11.49) 
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The suifajce foices and body forces are combined to generate the total force equations; 



(11.52) 


We now have the necessary expressions to set np the momentum equations for the 
geophysical fluid flow model. 


5, The Momentum Equation for Geophysical Fluid 
Flow 

Applying Newton’s Second Law (Fi = moi = pOidV) and using (11.45) and 
(11.52), we now write the momentum equations for the geophysical flow model: 

/Die \ dp dr^ dr^y dr^^ 

pa^ — p — f'uj — —— - 1 - - 1 - - 1 - 


dz dz 


By 


dz ’ 


pOy = P 




dr^ dryy dry^ 

dy dz dy dz ^ 


(11.53) 


pa 2 — 


0 




dr., 


dr,,. dr,. 


— — P9-\- ~ + + 


dz ' dy ' dz 

These equations are simplified using the assumptions described earlier: 


• The Fluid is Incompressible-. Density is independent of pressure, and therefore 
the model is uncoupled from thermodynamic considerations. 

• The Fluid Inuiscid'. All viscous forces are equal to zero. 

• The Fluid is Homogeneous-. We need not deal with the complexities of density 
stratification (this simplification will be relaxed in Chapter V). 

• Centrifugal Forces are Negated by Gravity-. This allows simplihcation of the 
acceleration terms. 
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m L r \ This allows us to neglect teims that arise in the acceleration compo¬ 
nents that result from the curvature of the Earth (e.g, see Cushman-Rjoisin 
(1994) [Rjef, 28]), 

• The Flow Primarily Horizontal'. Velocity and acceleration terms in the ^- 
directdon are negligible, 


Since viscosity terms are negligible (11,53) becomes: 



(11,54) 


n 


D 


Expanding the operator — introduced in (II. 11) produces the momentum equations 
in simplified form; 


. 1 

x-momentum: — -\-w— P - fr= - - 7^, 

at az ay paz 


^-momentum: 


^-momentum: 


dv dv dv 1 dp 

37 + "^3“ + ^'37 + = —3-j 

at az ay pay 


n 

0 = — -^- 3 - 

pdj; 


(IL55) 


Since geophysical flow is “primarily horisontal^, the term w— does not appear in 

az 

the z- and ^-momentum equations. Equation (11,55) together with the continuity 
equation (11,6) are used to construct a shallow water model. 


6, Governing Equations for the Shallow Water Model 

The shallow water model is depicted in Figure 5, We dehne a yariable h as 
the height of the fluid above a reference level ^ = 0, The variable k varies with z, 
y, and i and describes the fluid action at the surface. The “bottom^ of the shallow 
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Figure 5, The Shallow Water Model 


water environment consists of a rigid surface whose heigjit that does not vary with 
time and is given hy ^ = ks{z,y)- W'e let J/ = k{x,y,i) — ks(z,y) where H is the 
depth of the fluid layer with respect to the bottom contour. If L is the horizontal 
scale, then by the shallow water assumption it is true that H L. 

It is now possible manipulate the governing equations for the Earth model to 
produce the shallow water model, Integrating the .z-component of the momentum 
equation (11.55) yields: 


p{x, y, z, t) = -pgz + p(x, y, t) 


(IL5S) 


On the surface, .z = pressure p equals some ambient pressure Pq. Therefore: 


p(r, y, z, t) = pg[h{x, y, t) - r] + Po- 


(11.57) 


Dropping the variable dependencies for brevily, (11.57) is rewritten as: 


p = pg{h - .r) + Po. 


(II.5S) 


From (11.58) we obtain pressure gradients in the horizontal z and y directions: 


dp dk dp 


dk 


az az ay cfy 


(11.59) 
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Thus we lewiite the x- and ^-momentum expiessions (11.55) as; 


(11.60) 


Equation 11.59 aJso inipilies that the horizontal acceleration components, and 
Oy, given by (11.45) and (11.55) are independent of z. It follows that the horizontal 
velocity components, v. and are independent of z. This allows us to uncouple (II.6) 
and solve for w{x, y, z, t); 


Su du du dk 

x-momentum; — + if— - fv = 

cft ax ay ax 


dv dv dv . dk 

^-momentum; — + = -g^r^ 

at ax ay ch/ 


w{x,y,z, t) 


dx dy 


+ w{x,y,t) 


(11.61) 


Now consider the flow along the bottom contour ksi^, y)- Since the contour is rigid, 
there is no normal flow. Any velocity in the z-direction is due to fluid flowing tan¬ 
gent to the bottom contour. Figure 6 depicts this consideration for the x— and z— 
directions. If is the velocity component in the z-directi on flowing tangent to the 



Figure 6. Flow Along the Bottom Contour in the n-Direction 


bottom contour, its contribution to is; 

V. tan(Q!) = V. 


dks 

dx ’ 


(11.62) 
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wheie a is the ajigle of the line tangent to ks aJid x-axis. SaniilaTly in the ^-direction 


the contribiitaon to is: 


r tan(/3) = r 


dks 


(11.63) 


wheie i3 is the angle of the hne tangent to and ^-axds. It follows then that; 


dks , dks 


w{z,y, z = ks,i) = + V 


(11.64) 


dx By 

Theiefoie from (11.61) and (11.64) we generate an expression for w along the bottom 


(11.65) 


( 11 . 66 ) 


contour: 

. ^ ■. (Bu (9i' \ BKb BKb 

A corresponding condition at the surface, z = k yields; 

Bk Bh Bk 

. Bk , . . 

Since the surface is not rigid, we pick: up — in the expression. Therefore, (^11.65) and 

Bt 

(11.66) generate; 

Bk Bk Bk ^ . fBu Bv\ Bks Bk^ . . 


which is simplified to; 

( 11 . 68 ) 

Equations (11.60) and (11.68) are the governing equations for the shallow water model. 
In its current form the model is non-hnear. We desire a linear form for subsequent 
investigations. 


d d 

^ ^ [i' ('j - M] = 0- 


7, Linearizing the Shallow Water Model 

To linearize the shallow water model we conduct a perturbation analysis on the 
I- and ^-momentum equations (11.60) and the vertical momentum equation (11.68). 
Perturbation analysis on the horizontal flow is accomplished by assuming that the 
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Figure 7. Shallow Wa,tei Model Cioss Sectdon with k = H hs 


I- ajid ^-components of velocity aje dominated by constant teims {U and V). Sn- 
peiimposed on these are small variations, t) and t). Mathematically 

stated: 

y. = ?7 -1- y*, where U is constant and y*" 0{5) U, 
y = V -\-y*, where V is constant and y* 0{S) V. 

Applying these to (11.60) yields: 


By* y.By* 


— + ^U + u)^ + ^V+^)^ + f^V + u) = -g^. 

Ignoring teiins of 0(^^) results in the following simplification; 

oh;* ^,olt;* ..oh;* ,, dh 


dv*'- ^.dv* dh 


(11.69) 


(11.70) 


(11.71) 


A similar perturbation approach is applied to the vertical momentum equation 
of the shallow-water model (II.6S). We let k(x, y, t) = hs(x, y)-\-H(x, y, t) as depicted 
in Figure 7, Applying this to (II.6S) yields: 
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(11.72) 


&t 




% 


(yH) = 0 . 


wlieie H = k — ks- We lefiiie H{z,y,i) fuithei into two parts as depicted in Figure S. 
The component, Ho{z,y) is resultant from the constant velocity terms U and V 



Figure S. Shallow Water Model Cross Section with k = Hq -\- 7] -\- 


introduced in (11.69). If U and V are zero then Ho is constant. Likewise, if the 
bottom is fiat, then there is no velocity' in the vertical direction along the bottom 
that results from a non-zero U and V, and again Ho is constant. Superimposed on 
Ho is a small amplitudinal variation that represents the wave action and is given by 
7](jc,y,t). Mathematically stated this is; 

H(j:, y, t) = J/o(x, y) -\- 7j(z, y, t) where tj ^ 0(6) < Ho (11.73) 

Applying these terms, as well as the perturbation terms introduced for the horizontal 
velocity allows us to rewrite (11.72) as follows; 

^ + ^ [([/ + + >7)] + ^ l(v + + rj)] = 0. (IL74) 

Neglecting terms of gives us; 

^ ^ + v)+ ^ + i?) + = 0. (11.75) 
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Since U ajid V aje constant, 11.75 "becomes; 


I + + ,) + + ,) + = 0. (11.76) 


This can "be fnithei simphhed if the "bottom is flat (e.g, Ho is constant); 


(IL77) 


Finally, if we assume that V and U are zero (e.g. no advection), the equation becomes; 


dn „ /Sir' Si'’\ „ 


(11.73) 


This is the simplest form of the line arize d vertical motion component of the shallow 
water model. 

We now revisit the horizontal flow equations (11.71). By including the refine¬ 
ments and perturbations on k, (11.71) becomes; 


Su* T-^* r-T- f ^Ho ^r^ \ 




(11.79) 


Since we have assumed no advection (U,V = 0 and therefore Ho is constant) and that 
the bottom is flat (ks = 0), the linearized horizontal flow component of the shallow 
water equation can be stated as: 

(11.80) 

Equations 11.78 and 11.80 are the governing equations for the linearized shallow water 
model, They will be further refined in the following section, 
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8, Deriving the Klein-Gordon Equation 

We now c01111x116 lineaiized components of the shallow watei model to obtain 
the well known Klein-Gordon, oi dispersive wave equation. A step-by-step derivation 
follows. 

Step 1 : Perform the following substitutions to (II.73) and (11.30): 


ii = and r = 


The result is: 


^ - /r = 

Bv TV ^ 

Bn dii di 


(11,81) 

(ILS2) 

(11,83) 


Step 2 : Assume / is constant and take the partial derivative of (11.31) with respect 
to x\ 

B (Bu\ ^Bi- ^ 

and the partial derivative of (11.32) with respect to y\ 

a(ai\_,.aii_ ^a^ij 
ft/la J ft/ ^ °ay^’ 


and add the resulting equations: 

B f Bit f Bi 

Step 3 : Take the partial derivative of (11.31) with respect to y\ 

B fBii\ Bi B fBr]\ 

By ^ By ’ 

and the partial derivative of (11.32) with respect to x\ 

^ \ j. ^ \ 


(11,84) 
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(11.85) 


and subtract th 0 I 0 sulting 0 q'uatdons; 

d /du (9i'\ f (^ 


Step 4 : Tals t}i 0 partial d 0 iivativ 0 of (11.84) witli T 0 sp 0 ct to t and T 0 aTiang 0 tli 0 
I 0 sulting t 0 ims; 


/ 


Bt 


di; 

dz 


9y. 


I cm \ cf fr. \ 




Bz 


Bv 


s 


Mnltipily (11.85) by / and T 0 aTiang 0 tli 0 T 0 sulting t 0 Tins: 


-/ 


B (Bv 
(9t 


% 


= f 


Bv (9ii\ 
By Bz) ' 


Summing tli0S0 yi0lds: 



Bu Bv 
Bz By 




= 0 . 


Step 5: Using (11.83) i 0 wiit 0 (11.86) as: 


Bt 


B^rj 

Bt^ 


+ frj-gHo^^Tj 


= 0 . 


(II.S6) 


(11.87) 


Step 0: Int 0 gTating (11.87) with I 0 sp 0 ct to t and l 0 tting Cq = gHo yi 0 lds: 

( 11 . 88 ) 

wh 0 i 0 S{z, ^) is an arbritrary function of int 0 rgration. Equation 11.88, th 0 lin 0 aT 
inhomog 0 n 0 O'us form of th 0 Kbin-Goidon 0 qiiation, is a T 0 stat 0 m 0 nt of th 0 Iin 0 arb 0 d 
shallow wat 0 i 0 q'uation. It also d 0 SCTib 0 s oth 0 i b 0 haviors such as lateral vibrations 
of membrane strips and acoustic pressure waves in dispersive media [Ref. 24]. We 
continue with analytical considerations of the Klein-Gordon equation and the concept 
of dispersive waves in the next section. 


Bf^ 


+A=5(t,^), 


9. Analytic Considerations of the Klein-Gordon 
Equation 

Consider a homogeneous form of the Klein-Gordon Equation in one dimension; 


Sr2 


+ f^V = 0 - 


(11.89) 
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This equation has the folio wing solution; 


77(1, t) = exp[i(fci — wi)], 


(11,90) 


wheie: 

u'‘ = C^k^ + f\ (IL91) 

Here ti; is the angular frequency and k is the wave nunihei. Solving for k one obtains 
the dispersion equation for the Klein- Got don equation; 


k = 


- P 

Co 


(11.92) 


The phase velocity of the wave (e.g. the speed of a wave crest) is given by; 

= I- 


(IL93) 


Group velocity is the velocity at which the wave energy propagates and is given by; 


~ dk 


(11.94) 


If the phase and group velocities are not equal, the wave is dispersive and the wave 
shape deforms as it travels. Using (11.91) we have; 


k 


’t'-n — 




kCl 


(11,95) 


If / ^ 0, then I'p ^ Vg, and the solution to the Klein-Gordon equation produces 
dispersive waves. If / = 0, then ip = Vg = Co and the resulting wave is non- 
dispersdve. With respect to the Earth model, / is the Coriolis parameter given by 
(11.43). The magnitude of / increases as we go north or south from the equator. Thus 
dispersion effects will increase away from the equator. Hence the rotating Earth, with 
the exception of the equator where / = 0, is a dispersive environment with respect 
to the shallow-water model. 

This form of the Klein-Gordon equation is used in initial investigations of the 
Higdon NRHC because it is a relatively simple mathematical model of dispersive wave 
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"behavioi. These effects aje an uitegiail pajt in the development of Higdon NRBC’s. 
The simplicity of the equation also makes possible comparisons between numerical and 
analytical solutions and provides for easy testing of proposed boundary conditions. 
Chapter 2 follows with a description of the Higdon NRBC. 
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III. HIGDON NON-REFLECTING 
BOUNDARY CONDITIONS 

Oui goal is to accurately desciibe the piopagatdon of dispeisive waves in a 
“manage aid e^ domain. Snch a domain is one that will not over whelm computei 
capabilities when applying mimeiical techniques to its inteiior and boundaries. The 
actual domain in which the wave travels is much larger, or perhaps infinite. To restrict 
ourselves to a smaller domain of interest, artificial boundaries are constructed that 
allow waves that impinge upon them to freely pass. 

Constructing the mathematical analog for a NRBC, is elusive. Schemes exist 
that allow for the total absorption of non-dispersive waves. However, when dealing 
with dispersive waves, the results are less than perfect, and spurious reflections occur 
at the artificial boundary. Hence the NRBC problem is one of optimization (e.g. min- 
imisng unwanted reflection thus allowing most of the wave’s energy to pass). The 
Higdon NRBC is the focus for our investigation because it exhibits numerous advan¬ 
tages. We initially use a model that requires a single artificial boundary. From this 
we develop the details of the Higdon NRBC. 

A. THE SEMI-INFINITE WAVE GUIDE PROBLEM 

The semi-infinite wave guide problem provides a vehicle to explore the proper¬ 
ties of the Higdon NRBC and is depicted in Figure 9. A Cartesian coordinate system 

r 


Figure 9. Semi-Infinite Wave Guide 
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(x,^) is introd'uc0d such that the wave-gnide is parallel to the x-direction. The width 
of the wave-guide is denoted hv In the wave-guide we consider the inhomogeneous 
Klein-Gordon equation (II.8S): 


5=77 




5 , 


where, 17 is the unknown wave held, iGq is the given reference wave speed, / isthe given 
dispersion parameter, and 5 is a given wave source function. Cq and / are functions 
of location, but it is assumed that outside a hnite r^on they do not depend on x. 
The wave source S is a function of location and time, but is assumed to have local 
support. 

Dirichlet or Neumann boundary conditions are specihed on the north and 
south boundaries, Fs aJid Fjv: 


77 = 0 


or 



(111.1) 


In acoustics, these correspond to “soft waH^ and “hard waH^ conditions, respectively. 
A Dirichlet boundary condition is prescribed on the west boundary 


= VKv{y,-i^), (III. 2 ) 

where 7]w(y,i) is an incoming wave function. As x 00, the solution is bounded 
and does not include incoming waves. To complete the problem statement, the initial 
conditions; 

rj{x,y,[i) = T^, ^{x,y, 0 ) = Wi:i, (III. 3 ) 

are given at time t = 0 for the entire domain. We assume that the functions 170 and 
Wo have local support. 

We now truncate the semi-inhnite domain by introducing an artihcial ea^ 
boundary S at t which we call F^. This boundary divides the original semi- 

inhnite domain into two subdomains: an exterior domain I?, and a hnite computa¬ 
tional domain Fl bounded by Fjv, F5, F^;, and F^.-'. We chose the location of Fj; 
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such that the entire support of 5, t/oj ^ is in Thus on Tj; the homogeneous 
counterpart of the Klein-Gordon equation holds: 




-C^V^77 + /^77 = 0. 


On Tn, Ci and are ^-dependent (or, as a special case constant). 


Ba THE HIGDON NON-RE ELECTING BOUNDARY 
CONDITION 

The Higdon condition for non-dispersive acoustic and elastic waves was pre¬ 
sented and analyzed in a sequence of papers (see e.g. [Rjef. 17] - [Ref. 20]). These were 
later extended to include dispersive waves [Ref. 21]. To obtain a well-posed problem 
for the finite domain (see Figure 9) we impose a reformulation of the Higdon NRBC 
[Ref. 21] on to reduce spurious wave refection. 

The Higdon NRBC is obtained hy composing simple first-order differential 
operators. For example, a Higdon NRBC of order J (denoted hy Hj) is: 


(III .4) 


In this section we will show that Higdon NRRC’s have several advantages including: 

• Their reflection coefficients are easily determined. 

• They are exact for all waves that propagate in an j-direction with phase speeds 
equal to either of Cl through Cj. 

• They constitute a sequence of conditions of increasing order. No asymptotic 
approximation is involved in their construction, enabling one to obtain solu¬ 
tions with unlimited accuracy. 

• They are robust. Reflection coefficients become smaller as the order J in¬ 
creases. A good choice of C^’s leads to better accuracy for a smaller J, but 
reductions in spurious reflection can still be obtained with non-optimal C/s 
hy simply increasing J. 

• They are very general applying to a variety of wave problems in one, two, 
or three-dimensional configurations. Moreover, they can be used for wave 
problems in dispersive and stratified media. 
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To discav 0 T the motivatioii foi the foim introduced above (III.4), consider a 
possihle solution to the semi-infinite wave guide program suggested hy ( 11 . 90 ): 

17 = AYr, (y) exp [ifc( j - C^t -\- c)], (HI .5) 

where A is the wave amplitude, is the phase, is the t- component wave number, 
w is the wave frequency, and C^ = — is the phase velocity, Yniv) is determined from 
the dependency of Co and / on Rjecah from the Klein Gordon equation derivation 
that Co = where Hoix^y) is described in Figure S. Since this analysis assumes 
a flat bottom, J/o, and therefore Co, is constant. Also recall that /, the dispersion 
or Coriolis parameter, varied with the latitude ^ (11.43). With respect to the Earth 
model (Figure 2), 0 was determined to be a function of y only. Therefore, in the 
shallow water model, / is a function of y only. Now consider the real part of (111,5); 


Tj = AYr^{y) cos[fc(j - C^t -\- <■)]. 


(Ill. 6 ) 


Substitution into (III.4) generates: 


-4l^(^)fcsin [^(t — T V')] 


j 





(III.7) 


which implies: 

n(c, - c^)= 0. (iii.s) 

3=1 

Thus the Higdon NRHC is exact (e.g. no portion of the wave is reflected) at F^: if the 
phase speed C^ matches any of the chosen Higdon parameters Cj. 

Typically solutions to the Klein-Gordon equation consist of an infinite number 
of waves, whereas, we are limited to the selection of a finite number of parameters. 
Therefore, in most cases C^ ^ Cj. We can still validate (III.S) by assuming that the 
impinging wave splits at F^^ into a reflected and passing wave. The magnitude of 
the reflected wave is easily analyzed. Consider a simplified form of (III. 6 ) in which 
Yr,{y) = 1 and C = 0 : 

7] = A cos[fc(i — Ca;i)], (III,9) 
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and a fiist-ordei Higdon NRBC; 


_5 

dt 



’ 7 = 0 , 


(III.10) 


wheie ^ Cj. The oiiginal wave impinging on is the sum of the wave that 
passes thiough F^; and the wave that lehects hack into the domain. Mathematically 
stated: 


^ COs[fc( J — COs[ —-1- Ca-t)] -1- COs[fc(j — Ca-t)]. (Ill,11) 

wheie and Ap are the amplitudes of the lehected and passing waves lespectively. 
Note that lehection affects the wave hy leversing its direction of travel, or mathe¬ 
matically hy reversing the sign of the wave number k. The wave frequency ti; for 
the reflected wave remains unchanged. If any reflection occurs, then 0 < |^f.| 1-41. 

With r^ards to the passing wave, the wave number k and phase speed remain 
unchanged. However, |^| will be reduced with respect to \A\. Substituting the 
right-hand side of (III. 11) into (III. 10) yields: 


^,(C, + C^)-^p(C,-C^) = 0. 


( 111 , 12 ) 


Using this equation we define the reflection coefficient R for a first-order Higdon 


NREC to be: 


This yields: 




- a 







(III.13) 


(III.14) 


We see from (III. 13) and (III.14) that |^t,| < |^| when ^ 0. Note that R —1 as 


Ca; —0 implying that waves with low phase speeds will result in maximum reflection 
of \Ar\ = \Ap\ = .51^1 at the artificial boundary. This circumstance is mitigated in 
that these same low-speed waves reach the boundary at times that are often outside 
the scope of the problem. 
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The lehection coefficient for a J^^-oidei Higdon NHBC is easily derived using 
the above techniques: 

j 

B^n 

j=i 

This expression represents the product of J factors that are less than one. Hence sim¬ 
ply increasing the order J reduces the amplitude of the reflecting wave. Theoretically, 
one could reduce the amplitude of the reflected wave at to zero without r^ard for 
the value of Cj by letting J oo. Practically, however, we are limited by computer 
capability. We can only select a finite number of C/s and therefore must tolerate 
some spurious reflection. Fortunately, can by reduced significantly by intelligently 
selecting C/s. Strategies for this are discussed in the next section. 


Cj + C* 


(III.15) 


C. DETERMINING OPTIMAL VALUES FOR THE 
HIGDON PARAMETERS 

Existing literature offers no analytical means to optimize the value of the Hig¬ 
don parameters Oj for a dispersive wave. Three general methods have been suggested: 
(1) a-priori selection, (2) computer automated selection, and (3) dynamic selection. 
The first two methods are explored in this section. 

The first general method was suggested by Higdon (see e.g. [Rjef. 17] - [Hef. 
21]) and selects C/s using an “educated guess^, One takes advantage of the properties 
of the reflection coefficient J? and parameters of the Klein-Gordon equation to accom¬ 
plish this. A second general method chooses Cy’s automatically by computer code 
as a preprocess. These methods typically use information about the interior wave to 
select Oj’s that minimize reflection. Givoh and Neta suggest a simple approach in 
which Cj’s are determined from wave numbers that are evenly distributed over the 
span of maximum resolvable wave numbers. In another approach they recommend 
using the minimax method to pick k values [Rjef. 24]. To test these methods, the 
“oscilloscope^ method is developed to fine tune the C/s and ininiinize the reflection 
of a known wave. Theoretically, this procedure produces the best result, but it is too 


34 



0xpGnsdvG computationally to sgtvg as a prGpiocGss. It is, Iiowgvgi, a nsGlul mGasuiG 
of tllG GflGCtivGILGSS of OtllGT SUggGStGtl mGtIlods. 

To quantitativGly compaiG tliG diflGiGiit scliGmGS, SGVGial mimGiical GxpGii- 
mGnts ajG con ductGd foi a SGmi-in£nitG onG-dimGnsional domain [0, oo) on tliG T-aods. 
This domain allows wavGS to travGl in thG x-dir action only, Ghminating a need to con- 
sidai gGoniGtric dispaision which occuis in the sami-infinitG channal, It also allows us 
to datGimuiG analytically the affacts of Higdon paramatais Cj on tha i ah action coaffi- 
cient Tha wavas in tha domain ara go vain ad by tha ona-dimansion al homoganaous 
Klain-Goidon aquation: 


, ^2 


+ A = 0. 


(III.16) 


An artificial boundary is placed at 4 tt. A continuous wave is assumed to exist 
inside the truncated domain with the following initial conditions: 


75 


77(x, 0 ) = 51 ! F? ~ 


fc=L 


(111,17) 


The solution is assumed to be of the form: 

75 

i 7 (x, t) = ^ Afc cos(fcx — ujt), (111,18) 

fc=L 


where the dispersion relation: 


(111,19) 


is necessary to satisfy (111,16), Substituting (111,19) into (III.18) yields: 

7]{x, = T . (111,20) 

The resulting dispersive wave is somewhat contrived since the number of wave num¬ 
bers k is finite. However, it serves for the current purpose to analyze various schemes 
to select Higdon parameters, Rjewriting (111,20) generates: 

= I] 4(111,21) 
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wheie; 


^ + r 

k k 


(IIL22) 


Note that Ck ^ Co, a fact that is important later when estimating C^’s. Equation 
(III.21) is used to determine the reflection coefficient R\ 


„ Ar /r Cj-a A 


^=L C'j + Cfc ^ Ic 2J,2 ^ p 


(III.23) 


Thus the equation for the reflected wave is; 


£ h (m) ^). 


(III.24) 


which can be use to determine the error ||e'(t)|| resulting from a specified Higdon 
NHBC. Note that, by (III.23), /i is a function of k and must appear inside the 
summation. We calculate the error by taking the norm of the reflected wave from 0 
to r^;. Since the experimental data is discrete, the equation for the 2-norm of the 


reflected wave at time t is: 


k(i)ll = 




(IIL25) 


where N^, the number of elements in the x-vector, is determined by the refinement 
of the grid. In each experiment we set Co = 2 and / = 1. The error measurement is 
taken at time t = 100 time units. 


1. Experiment One: 

In the first experiment, we disregard the wave generated inside the domain 
and offer orur best guess for determining C^’s. As shown by (III.22), > Co- Fur¬ 

thermore, from (III.23), we see that for |Cofc| >> |/|, 


n (7 . Ca 


(IIL26) 


It is therefore logical to let Cj = Co for j = For a first order Higdon 

NRHC (Hi), this is equivalent to the Sommerfeld condition which is utilized widely 
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in meteoTology [Rjef. 24]. By simply incieasuig the oidei J of the Higdon NELBC, 
(III.15) TGVGaJ 0 d that thG iGflGction coGfliciGnt is TGdncGd. WliGn tvg iGt Co=2 and 
/=1, OUT GxamplG wavG (III.20) bGcomGs: 


75 


77( J, t) = 1] CQs(fcl - V4fc^ + 1 t). 

fc=L ^ 

The J*^-ordei leilectioii coefE.cieiit (111,23) foi this wave is: 


(III,27) 


TO-' 

j=l 


(111.28) 


2k + + 1 

Using this Gqnation wg gGnGiatG thG iGUGctGd wavG iGsuIting from thG boundary con¬ 
dition. FiguiG 10 (left) displays thG iGflGctGd wavG for a fiist-ordGi (Hi) throngji 
thiid-oidGi (J/j) Higdon NRBC. For through Hg thG amphtudG of thG rGflGctGd 
wavG has dGcrGasGd to thG point that it is not visiblG imlGss a smaJlGr scalG is nsGd. 
FigurG 10 (right) displays thG norms of thG rGflGctGd wavGS, which ajG talcGn to bG thG 
mGasuTGmGnt of Grror for thG booindary condition. This hgurG quantifies the rapid 
drop off in the refiected or error norm for increasing J. 



Figure 10. Left; Ekperiment la, at t=100; Cj = Cq with Cq=2 and /=! (Solid 
Line Depicts Hi. Dotted Line Depicts H 2 Plot). Right: ||i 7 h || at t=100 for Hi througji 
H^. 


We might refine this method even further by considering the following argu¬ 
ment; Although (III.23) reveals Cj ^ Co ^ k gets large, (III.28) suggests that R is 
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smaJl foT large k anyway. Therefore by selecting Cj = Co, we select Higdon coeffi¬ 
cients that damp out waves that had a low relection coefficient Ji, while waves with 
low wave numbers, and therefore higjiei reflection coefficients, are not considered in 
the estimation of Cj. Therefore, using (III.22), we consider: 


Cj = 


^Cq T f‘ 

kmiri 


(III.29) 


where, from onr experiment parameters, kmin = 1- Hoi this example Cj = 
Improved results are reported in Figure 11. However, one has presupposed some 
knowledge of the wave action inside the finite domain in order to estimate k^ir^. 



Figure 11. 


||77ji|| at t=100; 








V5 (Cc=2,/=1) 


2. Experiment Two: Cj Determined from Wave Num¬ 
bers Distributed Evenly over 

In the second experiment, we selected Cj’s by evenly distributing wave num¬ 
bers over the interval fcmoa;]. All other parameters used are the same as those 

used in the first experiment. Here, fcmmj "the minimum significant wave number, might 
be determined by some advanced knowledge of the internal wave, or by the time scale 
of the problem (e.g. kmiri describes a wave that moves so slowly that it never reaches 
the boundary). For this experiment, k^i^^ = 1. With regards to k^^, assuming 10 
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grid points pei wav 0 length aje necessary to resolve any given wave, a reasonable 
estiniate is: 


(IIL30) 

where Az is the grid spacing. The k values are then given by: 

kj, = ^ = (IIL31) 

where J is the order of the Higdon NRBC and kj. is the fc-valne for a J^^-order 
boundary condition. Equation (III.30) suggests that fcmaa: = 50 for this experiment. 
After determining the kjs, we now use (III.22) to calculate the Cj’s: 




yCpfej, H- 


(III.32) 


where Cj. is the Higdon parameter for a J^^-order boundary condition. The 
resulting C/s for this experiment are: 


Hi : 

2.2361 








H, : 

2.0001 

2.2361 







H^ ; 

2.0001 

2.0004 

2.2361 






H^- 

2.0001 

2.0002 

2 .000S 

2.2361 





H, : 

2.0001 

2.0002 

2.0004 

2.0014 

2.2361 




Fa : 

2.0001 

2.0002 

2.0003 

2.0006 

2.0021 

2.2361 



H, : 

2.0001 

2.0001 

2.0002 

2.0004 

2.0003 

2.0030 

2.2361 


Hs : 

2.0001 

2.0001 

2.0002 

2.0003 

2.0005 

2.0011 

2.0039 

2.2361 

Fg : 

2.0001 

2.0001 

2.0002 

2.0002 

2.0004 

2.0007 

2.0014 

2.0049 


The error measurements are reported in Figure 12 and represent a significant im¬ 
provement over those reported in the first experiment. We now seek to improve the 
estimation of C^’s using a method suggested by Neta and Givoh [Ref. 24]. 


3. Experiment Three: Estimating from Wave 

Numbers Obtained using the Minimax Formula 

The third experiment utilizes a method that computes C/s from wave num¬ 
bers obtained using the minimax formula based on the Chebyshev polynomial, an 
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Figure 12. Expeiiment 2, ||i7jj|| att=100: fioiu Iz Eivenly Distributed on 

(Co=2,/=l) 


approach adopted by Neta and Givoli [Rjef. 24]. We consider the nietbod for the one- 
dimensional problem presented by the current experiment in which wave numbers are 
given by: 






^ 2 


1 -|- COS 


( 2 ^- 1 ) 

2(J - 1) 


TT 




(III.33) 


where kj. is the wave number calculated for a order boundary condition. As 
before, the maximum resolvable wave number equals 50. The Higdon parameters 
are calculated using (III.32). The Higdon parameter Cj is set equal to Cq. The 
resulting Higdon parameters were; 


H, : 

2 



H, : 

2 

2.0002 


: 

2 

2.0001 

2.0007 


2 

2.0001 

2.0002 

H, ; 

2 

2.0001 

2.0001 

J/a ; 

2 

2.0001 

2.0001 

H7 : 

2 

2.0001 

2.0001 

Hs : 

2 

2.0001 

2.0001 

: 

2 

2.0001 

2.0001 


2.0015 



2.0003 

2.0026 


2.0002 

2.0005 

2.0041 

2.0002 

2.0003 

2.0007 

2.0001 

2.0002 

2.0004 

2.0001 

2.0002 

2.0002 
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2.0059 

2.0009 2.0030 
2.0004 2.0012 2.0104 














The eiioi noiins leported in Figure 13 (left) are not the predicted improvement over 
the results reported in the second experiment. However with a minor adjustment to 
the method, we can bring this about. Rather than setting Cj = Co, we use (III.29) 
and let: 



The Higdon parameters are now shifted slightly as follows: 


H, : 

2.2361 








H, : 

2.0002 

2.2361 







: 

2.0001 

2.0007 

2.2361 






: 

2.0001 

2.0002 

2.0015 

2.2361 





H, : 

2.0001 

2.0001 

2.0003 

2.0026 

2.2361 




Ho : 

2.0001 

2.0001 

2.0002 

2.0005 

2.0041 

2.2361 



H, : 

2.0001 

2.0001 

2.0002 

2.0003 

2.0007 

2.0059 

2.2361 


Hs : 

2.0001 

2.0001 

2.0001 

2.0002 

2.0004 

2.0009 

2.0030 

2.2361 

Hg : 

2.0001 

2.0001 

2.0001 

2.0002 

2.0002 

2.0004 

2.0012 

2.0104 2.2361 


The decrease in error measure reported in Figure 13 (right) indicates that though the 
change in the method is minor, the results are substantial and are a small improvement 
over the results of the second experiment. We now conduct a final experiment that 
will fine tune the results found in this experiment and ininimize the reflected wave. 

4. Experiment Four: A Procedure for Optimizing 

The fourth experiment utilizes a procedure for optdmiang C/s that can be 
utilized if an exact solution is known. It is computationally intensive, but provides 
a quantitative comparison for the methods introduced thus far. Imagine the Higdon 
parameters Cj to be controlled by the dials on an oscilloscope. A boundary condition 
of order J is represented by J dials. The screen on of the oscilloscope represents 
the finite domain and displays the reflected wave. An experimenter adjusts the first 
dial up or down in order to reduce the size of the reflected wave. He continues his 
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Figure 13. Left: Expeiiment 3a, ||i7jj|| at t=100 and = Cq. Right: Experiment 
3b, ||;7 j II --I ^ 


at t=100 and Cj = 




In both Gases, is Computed from 


Wave Numbers Determined using the Minimajc Formula Based on the Chebyshev 
Polynomial 


adjustments to the remaining dials minimizing the reflected wave each time. After 
finishing he repeats the process as many times as necessary until he can no longer 
reduce the displayed wave. The resulting dial settings are the Cj’s that minimize 
the reflection at the artificial boundary for a particular wave action in the domain of 
interest. 

This procedure was simulated with MATLAB and using parameters deter¬ 
mined from the results in the third experiment as the initial C^’s. The resulting 
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Higdon pajameters were; 


H, : 

2.2229 




H, : 

2.04S1 

2.2359 



: 

2.0170 

2.0607 

2.2361 



2.0075 

2.0260 

2.0615 

2.2361 

H, ; 

2.003S 

2.0135 

2.0275 

2.0615 

J/a ; 

2.0031 

2.0116 

2.0206 

2.0273 

H 7 : 

2.0021 

2.0075 

2.0151 

2.0276 

Hs : 

2.0012 

2.0047 

2.0094 

2.0156 

H. : 

2.0010 

2.0037 

2.0075 

2.0106 


2.2361 




2.0616 

2.2361 



2.0603 

2.0615 

2.2361 


2.0276 

2.0615 

2.0615 

2.2361 

2.0156 

2.0276 

2.0615 

2.0616 2.2361 


The eiioi measurements leported in Figure 14 show a significant deciease when com¬ 
pared to those geneiated "by the modified minimax method, dearly the “oscilloscope^ 
piocedme was effective in substantially reducing the magnitude of the reflected wave. 
It is also interesting to note that the maximum Cj for J /3 through J/g is based on 
which is 1 in this experiment. This has the effect of allowing the wave with the 
smallest wave number to pass through the boundary unhindered. This is significant 
because this wave has largest reflection. We can obtain the remaining wave num¬ 
bers that pass without reflection at the boundary using the above C^’s. These are 
calculated as follows: ,_ 


- 




/■ 


cl - 


(IIL34) 
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Using this foimula, the lesulting fcj’s fiom the oscilloscope piocednie are; 


H, ; 

1.0308 








H, ; 

2.2662 

1.0004 







; 

3.8237 

2.0139 

1 







5.7711 

3.0908 

2.0008 

1 





H, ; 

8.0799 

4.2890 

3.3035 

2.0008 

1 




ffa ; 

9.0027 

4.6387 

3.4745 

2.9870 

1.9992 

1 



; 

11.011 

5.7748 

4.0590 

3.0010 

2.0132 

2.0008 

1 


; 

14.333 

7.3193 

5.1511 

4.0014 

2.9986 

2.0009 

2.0002 

1 

; 

15.918 

8.1628 

5.7489 

4.8542 

3.9942 

3.0011 

2.0004 

1.9992 


This infoimation by itself is uninteresting, but is useful when compared to the kjs 
obtained by the modified minimax method which weie: 


Hi : 

1 




: 

35.355 

1 



Hz : 

46.794 

19.134 

1 



48.296 

35.355 

12.941 

1 

i/5 : 

49.039 

41.574 

27.779 

9.7545 

i/e : 

49.384 

44.550 

35.355 

22.700 

i/7 : 

49.572 

46.194 

39.668 

30.438 


49.686 

47.194 

42.336 

35.355 

i/, : 

49.759 

47.847 

44.096 

38.651 


1 


7.8217 

1 



19.134 

6.5263 

1 


26.602 

16.514 

5.5982 

1 

31.720 

23.570 

14.5142 

4.9009 1 


These lesults indicate that: 

• It is moie effective to select C^’s that coiiespond to lowei wave numbeis. 

• The span of wave numbeis generated by the irdnimax method is much gieatei 
than the optimized span indicated by the “oscilloscope’’ pioceduie. 

• The wave numbeis inside the minimax span increase too quickly when com¬ 
pared to the optimal span. 
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Figure 14. Ebcpeiiment 4, ||i7i?|| att=100: Cj Optimized using OsciUcEcope Pioceduie 


In the next chapter we explore the properties of the Higdon NRBC with respect 


to the two-dimensional infinite wave guide. 
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IV. DISCRETIZING THE SEMI-INFINITE 

WAVE GUIDE PROBLEM WITH ARTIFICIAL 

BOUNDARIES 

In this section Tve discretize the Klein- Got don equation. We also derive two 
forms of the discrete Higdon NRBC using the two- and three-point hactward-difference 
equations for the first-order derivative. A scheme proposed hry Givoh and Neta [Ref. 
22] is used to simplify this step. We then compare the two Higdon forms, as well as 
previously discussed schemes for selecting Cj, via numerical example, We accomplish 
this by numerically solving the semi-infinite wave guide problem (Figure 9), A spatial 
grid for j and y and a temporal grid for i is set up and he 17 value at a point (jp, y^) 
at time i is denoted bv: 


7^^=rj{^p,yg,t„). (IV ,1) 

The 17 values of the wave guide’s interior points are determined using the discrete 
Klein-Gordon equation. We then obtain 17 values on the north and south boundaries. 
Neumann conditions are imposed here and discretization schemes are required. Fi¬ 
nally, the 17 values on the east boundary are determined using the discrete form of 
the Higdon NRBC, 


A. DISCRETIZING THE KLEIN-GORDON EQUATION 

The homogeneous linearized shallow water with zero mean or Klein-Gordon 
equation is given by (II.S 8 ): 


at^ 


- + frj = 0 . 


Higdon proved that in the context of the scalar Klein-Gordon equation, discrete Hig¬ 
don NRBC’s are stable if such a standard second-order central-difference scheme is 
used to discretize (V.19) [Ref. 21]. Thus we use the following approximation for the 
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S0coiid d0Tivative: 


^"(^ 0 ) = ^ -h)- 2F(io) + F(io + fe)] + 


(IV,2) 


w}i0i0 kis t}i0 grid sdz0 and 0 {h?) is an appioxiniatdon 01101 of □id0T h?. W0 appiox- 
iniat0 tli0 Kl0in-GoTdon 0q'uatdon by; 

- (^) (’7p-U - 


“ (^) + Vpj,+i) + - 0 - 


Solving foi giv0s; 


+ +’7p^+i) + [2 - (/Ai)“]’7p4 - n 


(IV,3) 


(IV .4) 


W 0 ■US 0 this 0 q'uation to appToxiniat 0 th 0 int 0 iioi points of th 0 domain. Th 0 appiood- 
matdon 01101 is of oTd 0 i 0( Ax^, A^^, At^). Not 0 that th 0 sch 0 m 0 is two I 0 V 0 I in tim 0 
and thus on 0 I 0 q'uii 0 s an additional starting val'U 0 . This can b 0 accomplish 0 d using 
th 0 Matsuno sch 01110 [R 0 f. 33]. 


B. DISCRETIZING THE NORTH AND SOUTH 
BOUNDARIES 

If th 0 conditions on the north and south boundaries contain derivatives, then 
they too must be discretized. For example, consider the Neumann type condition; 

5 = "’ 

The south boundary is discretized by using a three-point forward-difference formula 


foi F'(ro): 


^'(^ 0 ) = + 4F(io + h)- F(rc + 2h) + 0{h^). (IV,5) 

Zfi 
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Thus the south boundajy condition is approxiniated by; 


J ^ 1-3^5 + Ks+i - < 5 + 2 ] + 0( Aj,') = 0, (IV ,S) 

wheie S is the index niajldng the south boundary, is the size of the ^-grid, and 
0 { Ay'^) is an approxiniation error of order Ay"^. We now approximate 17 on the south 
boundary: 

’S.S = ^<3+1 - ^’S.S + 2 ■ (IV .7) 

Similarly on the north boundary, three-point backward-difference formula for F^(j:o), 
an approximation of order Ay'^ for 17 is: 

^.jv ~ ~ 2 ^p-Jv-2 j 

where fV is the index marking the north boundary. 


C, THE DISCRETE FORM OF THE HIGDON NRBC 

The order Higdon NRBC Hj (III.4) is the product of J operators: 




(IV.9) 


To discretize this equation, use the two-point backward-difference formula for the first 
derivative: 

F'iio) = ^ IF(xo) - F{io - k)] + 0{h), (IV,10) 

where 0(h) is an approximation error of order k. The discrete form of (IV. 9 ) becomes: 

^ ‘'7j(i,y,t)-T^{i,y,t- At)\ 


n 

j=l L 




At 


^ 77 (j’, y, t) -I 7 (j’ - Az,y,t) 


(IV.ll) 


= 0 . 


The approximation error for this equation is of order 0 (Az, At). 
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Two special opeiatois are defined to further facilitate the discussion. The 
temporal shift operator, denoted hy Sf, affects a grid point as follows: 


= - At). 


(rv.i 2 ) 


Similarly, the spatial shift operator in the t- direction, denoted by S , affects a grid 


point as follows: 


Vpj! = = vi^p - S',, ip)- 


(rv,i3) 


The identity operator, denoted by /, has no effect when operating on the grid point, 
e.g, = 77 ^. Using these operators, the discrete form of the Higdon NRBC is: 




(rv.i4) 


where the index E corresponds to the location of the Higdon NRBC at t By 

collecting the operator terms /, Sf, and S“, (rV,14) can be rewritten as: 


Making the following substitutions: 


At At 

fl-i — 1 Ci ~7 . dp — —1, and Pi — —. 

^ Ai ^ Ai 


(rv. 16 ) 


allows us to rewrite (IV. 15) as: 


j 

P (aji + + ejS“) 17 ^, = 0 . 


(IV,17) 


Now consider a second order Higdon NRBC H 2 , which by using (fV.lT ) and 


expanding can be written as: 

[aiasJ^ -\-a1d2ISf -\-a1e2TS~ -\-d]_a 2 lS^ -\-d]_d 2 Sf'^ 

-\-d1e2SfS- -\-e1a2IS- -\-e1d2SfS- -\-e1e2 = 0 . 

In this form, Hj is represented as the summation of 3^^ terms: 


(fV.lg) 




(rv.i9) 
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pattern that allows ns to qniddy generate the discrete form of Hj for any order J. 
With regards to the index m( 2 ), a digit valne of 0, 1, or 2 implies that the coefficients 








Tatle III. Constructing P-m ajid Am from the Base Three Index 0221(3) 


Digit Position 

1 

2 

3 

4 

Digit Value 

0 

2 

2 

1 

Use Digit Value to 

Select Constant 




^ j 

Use Digit Position to 

Set Constant’s Subscript 

ai 

^2 

ez 

di 

Use Digit Value to 

Select Operator 

I 


S~^ 

Sf^ 


aj, dy, or ej respectively are part of the product that mate up The position of 

the digit determines the index j for each term, where we consider the left most digit 
to be in position 1, Similarly, the digits 0, 1, and 2 imply the application of the /, 
S^, or S~ operator respectively in the construction of ■ FoJ example, has 
SI terms. The index for the 25*^ term has a base 3 equivalent equal to 0221(3) (jiote 
that J digits are always used for Hj). Table III summarizes the application of the 
base three index 0221(3) construct Am and Pm- Therefore: 

-^25 —^ ^022L(3j “ J 

and: 

We can further generalize the operator equation as follows: 

= vT-U , (IV.24) 

where a and b are the number of times that the digits 1 and 2 appear respectively in 

^^(3)- 

Thus far we have indicated that this scheme will generate 3^^ terms, but this 
number is reduced considerably when one combines terms containing the same opera¬ 
tor. For example, in Table II the operators ISf, IS~, and S~Sf appear twice. Wben 
combined, the number of terms for the discretized form of H 2 is reduced by three. In 
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gGiL 0 Tal, teims will contain opeia-tois of the form ^ wheia 0 < a -1- 6 < J. The 


nnnibei of integei conibina,tions of a and 6 that satisfy the in equality is: 

. (J + 2)(J + 1) 

^ J ~ n J 


(rv.25) 


and lepresents the nnnibei of diffeient op eiatoi possibilities for Hj. This is significant, 
because it greatly reduces storage requirements as well as processing time when using 
higher order Higdon NRBC’s. For example, Hg is fully discretized using 55 vice 19,6S3 
(3®) terms, 

D, IMPROVED FORM FOR THE DISCRETE HIGDON 
NRBC 

In the previous section, we formulated a discrete form of the Higdon NRBC 
with an approximation error of order 0{k). The approximation error for the dis¬ 
cretized Klein-Gordon equation and for the north/south boundaries were of order 
0{ h}). This suggests that we might improve our scheme for discretizing the Higdon 
NRBC by using the three-point backward-difference formula which is given by: 

= 4[3F(io) - 4F(ro - /i) + F(xc - 2h)] + 0{h\ (rV,26) 

Using this in operator notation, the improved discrete form of the Higdon NRBC is: 


Collecting the opeiatoi teims /, ajid S~y (IV.27) can be rewiitten as: 


(rv,27) 


Making the following substitutions: 

At 4 1 4 At 1 At . . 

^ A^’ = “3’ = 3’ = “3^'= 3^'^’ 

allows us to rewrite (IV.23) as: 


1 At 


[] J T bjSf -\- CjSf^ T djS- -\- = 0. 


(IV.30) 
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In the Gxpajid 0 d foini, Hj is lepiesented as the summation of 5^^ teinis: 

5-^-L 

= 0, 


(rv.3i) 


n^O 


wlieie Am is a product of a^, bj, Cj, dj, and/oi Samilarly Pm is made up of a 

combination of operators /, Sf, S~ and/or 5“^, By letting; 

J 

Ao = '[[ai and Z* = AmPr^ri^^ , (TV.32) 

fn=L 

and noting that ^ ^ 0 we rewrite (IV.Sl) as: 

vh = (IV.33) 

Using a procedure similar to that employed hy Givoh and Neta [Ref. 22], we can 
easily sift through the algebraic complexities of (IV, 33), 

Consider again the coefficient Am aJid the operator Pm- W'e rewrite the index 
m in base 5 which we will refer to as . With regards to this index, a digit value of 
0, 1, 2, 3, and 4 imphes that the coefficients < 1 ^, 6 ^, Cy, dy, or respectively are part 
of the product that mahes up ■ The position of the digit determines the index j 
for each term, where the left most digit is in position 1. Similarly, the digits 0, 1, 2, 
3, and 4 imply the application of the /, Sf, Sf"^, S~, or operator respectively in 
the construction of ■ h'or example, has 3125 terms. The index of the 1454^^ 
term has a base 5 equivalent of 21304(5). Therefore: 

and: 

= c 1^*2 (^3 (lie's ■ 

We have indicated that this scheme will generate 5^^ terms, but this number 
is reduced considerably by combining terms with the same operator. Terms will 
contain operators of the form IS~^Sf^ where 0 < a-|-i < 2J. The number of integer 


combinations of a and b must be less than: 


. (2J+2)(2J+1) 

/ ^ J (J 

j=L 


(rv.34) 
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We note that ceitaiii opeiatoi com'banatiQns are not possible. Foi example, consideT 
Hi. The opeiatoT Sf^S~ is impossible because all three operators are required to 
construct the .5^“^ term leaving no possibility' for it to coexist with the S~ term. 
Taking into consideration these operators, the number of terms needed to describe 




(2J+1)^ + 1 

2 


Hence Hg is fully discretized using ISl vice 1,953,125 (5®) terms. 


(rv,35) 


E, NUMERICAL EXAMPLE: EMPLOYING HIGDON NRBC’S 
FOR A SINGLE WAVE PULSE AT Fw IN A SEMI¬ 
INFINITE WAVE GUIDE 


We now return to the semi-infinite wave guide depicted in Figure 9 of Sec¬ 
tion A. We let the wave guide width 6=5 and depth 0 = .1. The medium is 
dispersive with / = .5 and a gravitation acceleration of ^ = 10 is used. The initial 
values are zero everywhere, The boundary function on the west boundary is; 




.005© cos 


^(y-vo) 




(IV.36) 


0 otherwise , 

where r and to are the wave pulse center, radius and time duration respectively, 
The parameter values are set at ^ = 2.5, r = 1.5, and to = 0-5. 

An artificial boundary B is imposed at x = 5, defining the computational 
domain £7 as a 5 x 5 square. A mesh of 20 x 20 is used in £7. The two-point backward 
difference method as described in Section IV. Equation IV.23 is used to estimate the 
Higdon boundary values. The extended domain V for the reference solution 7]^^ is 
a 15 X 5 rectangle with a 60 x 20 mesh. An artificial boundary is imposed on V at 

X = 15, Any spurious reflections resulting from this boundary will not affect 

25 . . . 

if the run time is less than —. This the time it takes for a wave to travel from 
X = 0 (r^;) to X = 15 (Fijf.- for V) plus the time it takes for any refiection to travel 
back to X = 5 (Fh.- for £7). On the north and south boundaries we impose — = 0. 
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This “hard wall^ condition, which causes waves to bounce off the north and scmth 
boundaries, will result in additional geometric dispersion. 

Two cases with artificial boundaries, 17 ^and 17 ^^ 3 , are computed and jux¬ 
taposed to 77 ^^. For an NRBC with J = 4 is constructed with parameters 

Cj = {Co, Co, Co, Co}. Here Co = = 1 is the minimum wave speed. An error 

measurement error ||e'(t)||n on the truncated domain ^2 is obtained by using; 




^=L j=L 


(rv.37) 


where and Ny are the number of grid-points in the z- and ^-directions as de¬ 
termined by the grid spacing. For 77 ^^j, an NRBC with J = 1 and Cj = {Co} is 
constructed and its numerical solution is compared to 77 ^^ to obtain a second er¬ 
ror measurement. We rate the effectiveness of NRBC’s by comparing the values of 
||e'(i)||n. Note that ||e'(i)||n ^ a function of time and will start to increase as the wave 
impinges on the artificial boundary. 

Figures 15 through 19 show the solutions for 77 ^^, 17 ^^! and 77 ^^^^ at times 
t=l, 4, 5, 3 and 10 seconds. The top-left and -right plots depict 77 ^^^ on the truncated 
domain Q and extended domain V respectively (note that, though the extended 
domain is continuous, it is separated in the hgure so that 77 ___^ may be better contrasted 
with 77 ^^! and i 7 ^„, 3 ). The middle- and bottom-left plots correspond to 77 ^,,! and 77^,^3 
respectively. Two graphs on the center- and bottom-right present ||e'(i)||n for 77 ^^^ 1 
and 77 ^^^^j. Note that “HNRBC-2DR-IS-IL-UO-VO-TOP appears in the caption. This 
shorthand will be used throughout the paper to identify to identify the figures. It is 
defined as follows: 


• HNRBC: Higdon non-reflecting boundary condition, 

• tiDR; 71-dimensional rectangular domain, 

• TiS: HNRBC applied to 71 sides of domain, 

• 71L: 71-1 ayer problem. 
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• UmpTi; j’-directdon advectdon U =-.n (note; m denotes a niimis sign and p a 
dednial point) 

• Vpn: ^-diiection advection V = .n. 

• Tn; solution at t =n. 

Theiefoie “HNRBC-2DR-IS-IL-UO-VO-TOP lefeis to the solution at i = 1 of a prob¬ 
lem in which a Higdon NRBC as applied to one side of a two-dimensional rectangular 
domain consisting on a single layer with no advection. We will introduce multi-layer 
stratification and advection models in later chapters. 
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Figure 15. HNRBC-2DR-lS-lL-U0-VCtT01 with Sample Pulse Boundary Condition 


At t = 1 (Figure 15) the wave packet is still close to Fiv and no spurious 
reflections have occured. The plots for 17 ^^, are identical and ||f(t)||j^ 

for and is essentially 0 . 

At t = 4 (Figure 16) the leading edge of the wave padcet reaches B. A slight 
spurious reflection is measured for 17 ^^,,! and 17 ^,, 3 , but overall the three solutions are 
still very similar. At i = 5 (Figure 17) the wave packet crosses B. The 17 ^^,,! and 
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Figure 16. HNELBC-2DR-1S-1L-U0-VC^T04 with Sample Pulse Boundary Condition 


plots aje still indistinguishable from 17 ^^. At t = S and 10 (Figures IS and 19), 
most of the wave padcet has left fl and is now visible in the eurtended domain V. The 
17 ^^^ solution exhibits wave traces similar to those in 77 ^^, whereas the 77^^3 solution 
does not, The difference in scale for ||e'(^)||n leveals an improvement of one order 
of magnitude for 77 ^^!. This is quantitatively significant, but the qualitative results 
are of greater note. For 77 ^the wave’s energy has passed through B relatively 
unimpeded. On the other hand, 77 ^^,_j reveals visible spurious reflections resulting in 
a backwards-moving wave that pollutes the computational domain £7. 

This example illustrates the robustness of the Higdon NRBC. The perfor¬ 
mance of the NRBC is improved by simply increasing the order J. We now use this 
example to measure any improvement in the solution using the two-point Higdon 
NEbBC approximation derived in Section IV.C versus the three-point Higdon NRBC 
approximation derived in Section IV.D. In both cases Cj = {Co, Co, Co, Co} where 
Co = 1. The results in Figure 20 report a solution improvement of an order of about 
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Figure 17. HNELBC-2DR-1S-1L-U0-VC^T05 with Sample Pulse Boundary Condition 


1/3 when three-point approximations are used. We also use the current example’s 
parameters to compare results when the are chosen as a preprocess. An algo¬ 
rithm using the symmetric minimax formula (based on the Ghebyshey polynomial) 
proposed by Sommeijer et al, [Rjef, 32], In Case 1, the Higdon NRHC has an order 
J = 4 with Cj = {Co, Co, Co, Co} while in Case 2, J = 4 and the C/s are chosen 
as a preprocess. In both cases the 2"^ order approximation for the Higdon NRBC is 
used. The expected improvement is not revealed by Figure 21. This discrepancy will 
be checlced in a follow-up example. 
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Figure IS. HNELBC-2DR-1S-1L-U0-VC^T0S with Sample Pulse Boundary Condition 



Figure 19. EINRBC-2DEI-1S-1L-U0-V[^T10 with Sample Pulse Boundary Condition 
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Figure 20, ||e'(t)||n Plots foi 2-Poiiit vs, 3-Poiiit Higdon NHBC Appiaximataons 



Figure 21, ||e'(i)||n Plots for Cj = {Cq, Cq, Cq, Cq} vs, Cj Selected using a the 

Symmetric Minimax Formula with 3-Point Higdon NRHC Approximations 
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F. NUMERICAL EXAMPLE: SEMI-INFINITE WAVE GUIDE 
WITH A CONTINUOUS WAVE INTRODUCED AT 

The set up foi this example is the same as the pievious example with the 
following exceptions: Instead of a single pulse, the boundary function on is con¬ 
tinuous, linear combination of three waves with the form: 

I] ^T«cos cos(fc^j’(rV.3g) 

m=l \ D / 

whose parameters are selected a priori to be: 

= . 01 , . 01 , . 01 : 

Um = 1 , 2 , 2 : 

= .31, 1.37, 1.63. 

To satisfy the Klein-Gordon equation, the dispersion relation: 

< = (IV.39) 

must be invoiced to determine 

Thus km, using the parameters given for (IV.33), is: 

k^= . 11 , . 22 , 1 . 00 . 

Two cases, and aje again juxtaposed to 17 ^^ in Figure 22. For 

an NRBC with J = 5 and parameters Cj = {Cq, Cq, Cq, Cq, Cq} is used. For 
J = 2 and Cj = {Cq, Cq} is used. An error measurement ||e'(i)||n given by (IV.37) is 
determined. As expected, the higher-order Higdon NRBC results in a smaller ||e'(t)||n 
and therefore generates less spurious reflection at B. 

We now examine previous suppositions considered in Section III.G concerning 
the selection of the Higdon parameters C). 
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FiguTG 22. HNRBC-2DR-1S-1L-U0-V0-T10 with Contdimous Wave B 01111 dary Condi¬ 
tion 


• In Figuie 23 we compare Cj = {1,1, 1,1,1} vs. a 5^^-order Higdon NRBC 
where the are automatically selected using the minimax formula based 
on Chebyshev polynomials. The result favors the former, discounting previ¬ 
ous supposition that auto-selection of C^’s using this scheme might yield an 
improvement. 

• In Figure 24 we compare Cj = {1,1,1,1,1} to a 5^^-order NRBC where the 
tTj’s are evenly distributed from to -\- /^. Any improvement indicated 
by using the latter scheme is not significant. 

• In Figure 25 we compare Cj = (1,1,1,1,1} to a set of five C^’s that are evenly 
distributed from Co to the maximum resolvable wave number kma^ = —-—. 
Again the simpler scheme for choosing C/s is numerically superior. 

Based on results in this and the previous section we conclude that spurious reflections 
at the Higdon boundary are reduced by using: 


Cj’s equal to {Cq, ...Cq}. 
Higher order Higdon NRBC’s. 
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• thiee-point appioodmations foi the deiivatives in the Higdon NRBC. 





Figure 23, ||e'(t)||n for Cj = {1,1,1,1,1} vs, Cj Auto-Selected Using Mininiax For¬ 
mula Based on Ghebyshev Polynomials, 

We make one final comparison before closing. In Figure 26 a Higdon NRBC 
with J = 5 that uses a two-point approximation is juxtaposed to a Higdon NRBC with 
J = 4 that uses a three-point approximation. We see that the results are quite dose. 
However differences in computational effort is significant. From Section IV,C we know 
that to discretize a two-point Higdon NRBC with J = 5 required the generation of 
125 (3^) terms that were subsequently reduced to 21 terms when redundancies are 
combined (IV,25). From Section IV.D we determined that to discretize a three- 
point Higdon NRBC with J = 4 required the generation of 625 (5“^) terms that 
were subsequently reduced to 41 terms (IV.35). Since the discrete Higdon NRBC 
formula must be applied every time an artificial boundary grid point is encountered, 
it is apparent from this comparison that it is more efficient to use the two-point 
approximation with a higher Higdon order J. 
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FiguTG 24. ||e'(t)||n foi Cj = {1, 1 , 1 , 1 ,1 } vs. Cj Distri’b-uted Evenly fioni Co to 
+ p 



Figure 25. ||e'(t)||n foi Cj = {1,1,1,1,1} vs. Cj Distri'buted Evenly from Co to 

1 : — 
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Figure 26, ||e'(i)||n for Cj = {1,1,1,1,1} Usdng a Two-Point HNREC Approximation 
vs. Cj = (1,1,1,1} Using a Three-Point Higdon Approximation 
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V. 


A N-LAYER STRATIFIED DISPERSIVE 
WAVE MODEL 


We now lift the homogeneaus fluid assumption and develop equations to model 
geophysical flow in a stratified medium. A smtatle medium for a geophysical dynam¬ 
ics is the open ocean where fluid density p is affected by sahnity and temperature. 
Salinity changes are sfigfit in this environment, and temperature remains relatively 
Constantin the horizontal directions. However, temperature does change significantly 
in the vertical direction. Therefore we approximate pto be a function of s only. Since 
the fluid is assumed to be incompressible, p does not vary with pressure p. 

The linearized shallow water equations (11.73 and 11.30) were derived in part 
from the continuity equation for homogenous, incompressible fluids (II.6): 


du dv dm 


It was shown by (11.59) that v. and v are independent of s for a constant density 
fluid enabling us to uncouple m from v. and i' in (II.6), which after integrating yielded 
(II.Sl); 

(Sv{z,y,t) av{z,y,t)\ . 
w{x,y,z,t) = -z I- ^ -1 t). 

This critical step in the derivation is no longer possible when we assume that is 
dependent on z. 

We extricate ourselves from this conundrum by developing a layered shallow 
water approximation where p is constant in each layer (Fig. 27). Here it is assumed 
that the fluid is still incompressible and that density pi is constant in each layer Li, 
but varies in the different layers. In order for this stratification scheme to be stable, 
pi must be monotonically increasing downward [Rjef. 23]. Additionally we assume 
that there is no fluid mixing between layers. 
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Figure 27. A^-Layei ShaHow Water Model 


Refeiiiiig to the 7V-layei shallow water model, the pressure at any point in 
is determined from hydrostatic principles: 

^-L JV 

+ , (V.l) 

j=L j=i / 

where Po is a- constant ambient pressure at the surface ho and N is the total number 
of layers in the model. In (V.l), the first summation term is the contribrutdon to pi 
from the layers above Lj. The second summation term is the contribution to pi from 
the liquid column in Li. We use (II. SO) and (V.l) to obtain the horizontal momentum 
equations in L^: 

(V,2) 

where Tii, and Wi are the i-, y-, and z-components of velocily in Li. 




6 S 


Deiivation of the veitical momentum equation in is moie complex. Since 


pi is constant in we can uncouple the continuity equation foi 


,, (Sui{x,y,t) dri{z,y,t)\ , 

■mi{x,y,z,t) = -z [--1 ^V,3) 

wheie Wi is a veitical velocity component in Li. Foi "bievity we diop dependent 
vaiiahles fiom subsequent expressions. At the interface between Li-i and Li where 

JV 

"ths vertical speed component Wi (see similar discussion in Section II.B.6) 

j=i 

is; 

^ JV ^ JV o JV 

y—i j=i “ j=i 

This imphes that; 




&t ^ ^ dz 


(V,5) 


JV 


Similarly, at the interface between and Li^i where vertical speed 

J=-! + L 

component Wi is; 


^ JV ^ JV ^ JV 




(V,6) 


j=i+l 


which implies that; 


^jv d / ^ \ d / ^ \ 

* ■ «,£/'+s r i" 

Since Wi is independent of (V,5) and (V,7) must be equal, Therefore; 

(V.S) 

Equation (V.S) is the vertical momentum equation for Li. Together with (V.2) this 
completes the description of the fluid motion inside of the layer Z^. 

Before considering a numerical solution we linearize the governing equations 
for each layer. We assume that the and are dominated by constant terms 

Ui, V'i and ©^. Superimposed on these are small variations lij', r*, and of 0{6), i.e.; 

T Vi = Vi-[- V'i, and = ©^ T pi. (^ ■^) 


3 3 3 

-h + + -{x^h,) = 0, 
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Substntutirig these in (V.2) ajid (V.S) ajid neglecting teinis of 0(S'^) yields; 




Bt 


J=^ 


di 


I Tj I I f{Tj \ \ ^ ^-v in^ 


^ 4.V^' 4- +0 (^’^4- - 0 

If we assume that theie is no advection (e.g. = V] = 0), then (V.IO) lednces to; 



(V.ll) 


Equation (V.ll) is the linearized shallow water equation for the stratified iV-layer 
model with no advection. 


A. REDUCING THE STRATIFIED 7\r-LAYER MODEL 
TO THE KLEIN-GORDON FORM 

In order to combine linearized components of the stratified AT-layer model 
(V.ll) to a Klein-Gordon form, we use a step-by-step derivation that is similar to 
that used in Section II.B.S. 

Step 1 ; Perform the following substitutions in (V.ll); 

Hj = and 
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RecaJling that is independent of time, the lesult is: 


— ff — —-1- ^ 


j=t 


a />* S;, ^ , 

(5t7o <9uo 

— -I- —- A- —^ = 0. 
dt dx dy 


(V.12) 


(V.13) 


(V.14) 


Step 2 : Let / be constant. Talce the partial derivative of (V.12) with respect to x\ 

3r ( St j Si ' pi Si^ ^ Si^ J ’ 

and the partial derivative of (V.13) with respect to y\ 


;=! 


dy \dt j dy ^ V 


Add the resulting equations: 

d ^ dUi ^ dvi ^ ^ dvi dUi 

\_j = l r-L 

Step 3 : Take the partial derivative of (V.12) with respect to y\ 


St V Si + ^ 


(V.15) 


(9i/ (9t J dy ^ dy pi dx dx ^ 

and the partial derivative of (V.13) with respect to x\ 

(9t J dx ^dx pi dy ^=i ^ ^ 

Subtract the resulting equations: 

d ( dui \ ^ ( dvi dui \ 


(VIS) 


Step 4 : Take the partial derivative of (V.15) with respect to t and rearrange the 
resulting terms: 

ji.. 1 1 1 1 zZ „ "I" zZ ] ■ 


(9t V 5x dy J df^ \ dx dy 


.j=i pi 


J=^ 
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Multipily (V.16) hy f and learTangG the lesultdng teims: 


d /^ _ rlf 

dy ) dz ) 


Add the lesultdng equations; 


Step 5: Using (V,14) lewiite (V,17) as; 


(V,17) 


_a 


S"'?. , X2_ " 




+ f\ - 30iV’ 


Step 0; Integrating (V.IS) with respect to i yields; 


= 0 . 


(V.18) 



(V,19) 


where Si{z,y) is an arhitrary function that results from integration. Equation V.19, 
resembles the Klein-Gordon form reported in (II.S3), and is a restatement of the set 
of linearised (about a zero mean) shallow water equations (V.ll) generated for the 
AT-laver stratification model. 


B. DISCRETIZING THE KLEIN-GORDON FORM OF 
THE STRATIFIED A-LAYER MODEL 

A standard second-order central-difference scheme is used to discretize (V.19). 
On a spatial and temporal grid of orur choosing, we let 7]^^ be the FD approximation 
of r^i{z,y,i) at grid point (zp^y^j) and at time t,, in layer L^. Thus solving explicitly 
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foT yields: 



wlieie Coi = . We use tliis interior scheme in the numerical experiments pre¬ 

sented in the next section. The Higdon NRBC is applied independently to each layer 
and takes the form: 

21 ) 

The development of its discretization follows the discussion in Section IV.C. Hard wall 
or Neumann conditions are also applied independently to each layer and discretized 
accordingly. Since the discrete forms of the Higdon NHBC, Nuemarm condition, and 
the time stepping scheme (V.20) are exphcit, the whole scheme is explicit. 



C. NUMERICAL EXAMPLE: THE STRATIFIED N-LAYER 
MODEL 

Consider a geophysical process that occurs in a semi-infinite channel (Fig¬ 
ure 9). AH assumptions and simplifications used to derive the AT-layer model apply. 

A Cartesian coordinate system (x,^) is introduced such that the channel is parallel 
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to th& j’-diiGction. On the noitli ajid south honntlajiGS Fjv aJid wg spGcify thG 
NGumajin condition: 

^ = 0 for i = (V,22) 

dy 

On thG WGst "boundary F^J^/ wg piGsciilDG using a DirichlGt condition, i.G., 


= ^ = 1...N, (V.23) 


whGiG (y, t) is a givGn function for an incoming wavG or disturbance. At t — exj the 
solution is bounded and does not include any incoming Tvaves. The initial conditions 
are: 

0) = 0, ^ Q i=l...N. (V.24) 

To obtain a well-posed problem in a hnite domain we impose a high-order Higdon- 
NEFBC on the east boundary F^^. 

We now apply the new stratification scheme to a test problem using the semi- 

JV 

infinite wave-guide. A channel width 6=5 and depth d = ^ = .1 aje se- 

^=L 

lected (note that to satisfy the shallow water assumption, it must be true that depth 
d -^horizontal dimensions). The stratified medium is modeled with six layers. The 
layer thicknesses from top to bottom are 0^ = {.01, .01, .01, .01, .01, .05}. The density 
for each layer is given by p, = {1,1.05, 1.1,1.15,1.2,1.25}. A gravitational parameter 
^ = 10 and a dispersion parameter / = .5 is used. 

The boundary function 17 ^^. on F^y is stipulated to simulate two geophysical 
events. A surface disturbance, aldn to the wind acting on the ocean, is initiated in 
Li by setting: 




.0005 cos 


2.5) 


sin Ki 


0 < t < 12.5 , 


(V.25) 


\ 0 otherwise. 

Note that the maximum amplitude of the disturbance is small relative to the layer 
thickness ^]_ = .01 so that the validily of the model, which is based on perturbation 
analysis, is not violated. A second disturbance, simulating seismic activily on the 
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□ceaji flooi, is initiated in the Lq and given hy: 


I .001 

I 0 


if \y-2.b\ <1.5 L 6 < t < 7.5 , 
otherwise , 


(V.26) 


An other values for are zero, The simulation is run for 15 time units, 

The problem is solved for three different scenarios. First, an eurtended domain 
V is constructed using a 15 x 5 rectangle with a 60 x 20 mesh to compute a reference 
solution . Then two additional solutions, 17 ^and 17 ^^^ 3 , are computed on a 
truncated domain £7 in which an artificial boundary B imposed at x = 5, A 20 X 20 
mesh is used on the resulting 5x5 square so that the mesh for £7 and 17 are identical 
in the truncated region. For a Higdon NR.BC of order J = 5 with parameters 

Cj = {1, 1,1,1,1} is used for each layer. For J = 2 and Cj = {1,1} is used. 

Note that the value 17 in each case represents the total perturbation on the domain 
surface and is a combination of perturbation effects in each layer. 

The reference solution 17 ^^/ is then juxtaposed with and 17^^3 both graph¬ 
ically and quantitatively, The respective numerical solutions are used to produce filled 
contour plots that provide a “fingerprint^ for the resulting wave action, All contours 
and color schemes are relative to 17 ^,^^. These solutions are then used to obtain error 
measurements ||e'(t)||n which are calculated by: 


k(i)L = mi 

^=L j=L 




N^Ny 


(V,27) 


where and Ny are determined by grid spacing. In both cases, the error at the 
surface and for each layer interface is plotted versus time. Note that the size of 
17 precludes spurious reflections from polluting 17 ^/' Therefore ||e'(t)||^ serves as a 
measure of spurious reflection at B. 

At t = 7 (Figure 28), the surface disturbance in has populated £7 and 
is now visible in 17. In addition, the bottom disturbance has been initiated and is 
propagating in £7. Both l a^id 17^^,,3 exhibit wave traces similar to those in Tjref • 
However, the error measurement for 17 ^^! is an order of magnitude smaller than that 
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Figure 2S. HNREa2DR-lS-6L-U0-V0-T07 

□f 3 . This demonstrates that simply increasing the Higdon NRHC order J reduces 
spurious reflection. A description of the shorthand notation used in the figure caption 
was given in Section IV.E. 

At time t = 15 (Figure 29), the surface disturbance, which ended at t = 12.5, 
continues to propagate in and V. The bottom disturbance has successfully passed 
through without a significant increase in spurious reflection. The wave trace for 
17 ^^^ closely resembles that of rjref, however deviations iii t 7^^3 are now visible. This 
example demonstrates that a properly constructed Higdon NRHC can be used to 
restrict the domain of the AT-layer stratification model governed by the linearised 
shallow water equation. 
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Fig^e 29. HNREa2DR-lS-6L-U0-V0-T15 
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D. ADDITIONAL OBSERVATIONS FOR THE 
STRATIFIED N-LAYER MODEL 

Using th 0 same piolilem parameters as in the previous example, we now con¬ 
sider the behavior of the perturbation at each layer as predicted by the model. Filled 
contour plots are constructed to measure the magnitude of perturbation for the layers 
Li through Lq . All contours and color schemes are relative to the perturbation in . 
In Figure 30, we see that the action for both the surface and bottom initiated events 
has affected all layers, but most of the energy transmitted by the wave “sinks^ to the 
lower, denser layer Lq. This remains true at i = 15 (Figure 31) at which time both 
events have passed through the truncated domain fi. 



Figure 30. HNRBC-2DR-1S-6L-U0-V0-T07 Layer Perturbation Comparison 

We further consider model predictions for the behavior of the perturbation 
at each layer interface. Recall that the disturbance at a layer interface is the sum 
of the disturbances of all layers below the interface. Filled contour plots are now 
constructed for the surface and the five layer interfaces using data obtained from 
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FiguTG 31. HNRBC-2DR-1S-6L-U0-V0-T15 Layer Pertuibatdon Comparisoii 


AH contours and color sdienies are relative to the total surface perturbation. 
The first comparison is again at t = 7 (Figure 32). The contour plots reveal that the 
magnitude of the response to the surface event is damped with each successive layer 
interface downward. However, the wave appears to maintain its character at each layer 
interface as far as wave speed and geometric dispersion is concerned. We also note 
that the bottom layer event on the left side of the plot has immediately propagated 
to the surface. This is reasonable since the fluid is incompressible. However, as 
the wave propagates through Q, the effect on the lower layer interfaces once again 
dampens. At t = 15, the damping phenomenon in the lower layers is more pronounced 
(Figure 33). However, there is an exception on right side of the plot. In this regime, 
the perturbation is nearly zero and affected primarily by spurious reflection from B. 
Differences in each layer’s reflection may have caused this visible anomaly. 

For a final comparison, we consider a six-layer, two-layer, and single-layer 
model. For this we set up the parameters for two-layer model as follows: Si = 
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FiguTG 32. HNRBC-2DR-1S-6L-U0-VCI-T07 Layei Inteiface Peituibataon Compajison 


{.05, .05} and p. = {1.1,1.25}. Essentially we have conihined the nppei £ve layeis of 
the six-layer model and used the average density for the £rst layer of the trwo-layer 
model. The second layer of the trwo-layer model is the same as the bottom layer 
of the six-layer model. The parameters for the single-layer model are ^ = .1 and 
p = 1.175 (note that from (V.19) we conclude that density is not a parameter in the 
single-layer model and therefore the value for p is irrelevant ). The problem is again 
run for 15 time units, however the domain is extended to x = 15 for each model and a 
60 X 20 mesh is used. This ehminates artihcial boundary effects and comparisons can 
be made without concern for spurious rehection. All other parameters are the same. 
The results of the three models are presented for t = 15 (Figure 34). Comparing 
the contours of the models reveal that the surface effect is reduced as the number of 
layers increases. Hence, the single-layer representation tends to over-predict surface 
wave action if the medium is stratified. 

We now wish to apply the AT-layer model to an ocean environment. However, 
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FiguTG 33. HNRBC-2DR-1S-6L-U0-VCI-T15 Layei Inteiface Peituibataon Compajison 

befoie doing so, we must £ist extend the application of the Higdon NRBC to two 
dimensions. This will enable ns to model a “patch of ocean^ vice “a semi-infinite 
channel with hard walls^. 
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Figure 34, Model Comparison; HNaBa2Da-lS-6L-U0-V0-T15 \’s, HNaB02Da- 
ia2L-U0-V0-T15 ^’s. HNaBO2Da-iaiL-U0-V0-T15 
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VI. APPLYING HIGDON NRBC’S TO TWO 
OR MORE SIDES OF A TWO-DIMENSIONAL 
RECTANGULAR DOMAIN 


The equataons and concepts are now in place to construct finite domains in 
which Higdon NRBC’s are imposed at two or more boundaries. In this chapter 
we employ the use of several examples to check: the viabihty and the stability of 
schemes that utilize Higdon NRBC’s on two and four sides of the single- and multi¬ 
layer models, Special attention is given to the corners where two Higdon boundaries 
intersect, 


A, EXAMPLE ONE: A TWO-DIMENSIONAL SINGLE¬ 
LAYER SCHEME WITH HIGDON NRBC^S ON TWO 
SIDES 

For this example we employ a quarter plane shown in Figure 35, The semi- 

D 


Figure 35. The Semi-Infinite Quarter-Plane 


infinite domain V is bounded by F 5 and F^/ and is represented by a 10 x 10 square 
with a 40 X 40 mesh. The truncated domain is a 5 x 5 square with a 20 x 20 mesh 
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and bounded b^ Fjv, ^s, aJid Fiy. On both domains Az = .25, Ay = .25, and 

At = .025. A gravitational parametei ^ = 10 and a dispersion parameter / = .5 is 
used. A single-layer non-stratified medium with a thickness of d = 0.1 is considered. 
The problem is rnn for 15 time steps. 

The specification for the Higdon NEbBG on Fjv is given by: 




j=i 


^y. 


(VI.1) 


This equation is the ^-direction analog of (IV.9) which was used to describe the 
Higdon condition for A discussion equivalent to that in Section III.C applies in 
discretizing (VI.1). If Ax = Ay on the specified grid, and the same Cj’s are used for 
each boundary, then the parameters and operators used to discretize the boundary 
condition may also be used to discretize the Fjy boundary condition. To complete 
problem description, a Neumann condition is imposed at F 5 and given by; 



and a boundary function is imposed on Fi\y and is given by: 

' 5 

^ Z! 

fn=L 

0 

where the parameters and are selected a priori as 




n^Tr{y - 1.875) 
3.75 


sm 




a y < 3.75, 
otherwise , 
follows; 


(VI .2) 


j4„ = .001, .002, .001: 

TifY2 = 1, 3, 1; 

tiJm = -31, 1.37, 1.68, 


(VI .3) 


The restriction y < 3.75 is discussed in the next section. Homogeneous initial condi¬ 
tions apply. 

A solution and ^ computed for V and respectively. Filled contour 

plots are generated for each solution. In Figure 36, these solutions are reported at 
t = 2. The upper-right subplot displays the solution for on V. The upper-left 
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Figure 36. HNRBa2DR-2S-lL-U0-V0-T03; Wave Pulse Passes through Fjv 

subplot is a magniJicatioii of i 7 _^^ on the sub-domain £7 and represents a solution on £7 
as if no boundaries were present atj: = 5and^ = 5. The lower-left subplot reports 
the solution for oJi f7 where Higdon NRBC’s are imposed at Fjv and A 

qualitative comparison can be made between the two left-side subplots to determine 
the effectiveness of the NRBC in restricting the domain. A quantitative measurement 
of this comparison appears in the lower-right subplot, which reports an error measure 
||e(i)||n (See rV.37). Note that for computational purposes, Higdon NRBC’s with 
Cj = (1,1, 1 ,1 , 1 ,1 ,1} are imposed at ^ = 10 and t = 10 on 27. Any reflection from 
these extended boundaries of V will have httle or no effect on the domain of interest 
S7. 

In Figure 36 at f = 3, the wave pulse has reached Fjv. A visual inspection 
reveals no discernible spurious reflection. The magnitude of the error measurement 
is within acceptable norms (i.e. less than the error due to the discretization). 

Atf = 5 (Figure 37), the wave pulse has reached F^ and passes without visible 
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Figure 37. HNREa2DR-2S-lL-U0-V0-T05; Wave Pulse Passes tinougli Fe 


lefiectdon or undue increase in ||e'(t)||n. Note that a leading edge of the pulse is about 
to reach the intersection of Fjv and F^. At t = 15 (Figure 33), the wave has passed 
through the corner unimpeded with no adverse effect on ||e'(t)||n "ths left¬ 
side subplots are similar with no visible differences. The measured error ||e'(i)||n 
stabilized at approximately 2 x 10“'^. 

The corner, however, is still cause for concern, because two different ways to 
approximate the Higdon boundary grid points are possible. In the first method we 
initially evaluate the Fe gjid paints excluding the corner point. We then evaluate 
the Fjv grid points including the corner point. Hence the corner point is evaluated 
based on boundary values calculated for F^. This method was used in generating 
Figures 36 through 3^' I^ ^ second method, the grid points for Fjv (excluding the 
corner point) are fully evaluated £rst followed by the grid points of F^: (including the 
corner point). Now the corner point values are obtained from Fjv boundary values. 
R^ardless of the method used, the value of the corner point should be the same. This 
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FigiiTG 3S. HNREC-2DR-2S-1L-U0-V0-T15: Wave Pulse Passes tliiough Comer 


was tested and tlie results are compared via the error measurement pilot in Figure 39. 
The figure indicates that the two solutions are exactly the same. Hence no special 
considerations must he given when handling corner points that are the intersection 
of the Higdon boundaries. However, one must be careful to ensure that the corner 
point is not evaluated twice when designing algorithms to compute these boundary 
values. 
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Figure 39. Coinei Point Check:: Two Approaches to Evaluating HNEhBC Grid Points 
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B. INSTABILITIES IN THE 2-D, 2-SID ED HIGDON NRBC 
SCHEME 



Figure 40. HNRBO2DR-2S-1L-U0-V0-T03; Fjv InstaMty with Cj = (1,1,1,1,1} 


Despite the heartening results of the previous example, a stability problem 
does emerge along the entire northern NRBC if a non-zero boundary function 
Is extended to Fjv Recall the constraint y < 3.75 in (VI.4). This was specified in 
order to set up a buffer consisting of five zero valued grid points between the boundary 
function in the northern Higdon boundary. If this is not done, the boundary instability 
occurs and propagates through £7. For example, define the boundary function at Fw 
as; 

cos (VI.4) 

The grid point at ^ = 5 is immediately adjacent to Fjv. All other parameters with 
regards to the previous example remaining unchanged. The problem is run with J = 5 
and Cj = {1,1,1,1,1}. In Figure 40 at i = 3, we discern turbulence in the upper 
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left-hand cornei of £7 foi r^ca^i ■ Ihi Fignie 41 at t = 10, we see that the instability^ has 
piQpagated through Q. 
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Figure 41. HNRBa2DR-2S-lL-U0-V0-T10; Fjv Instability with Cj = {1,1,1,1,1} 


To explain the unstable behavior, recall that the 5^-order Higdon NRBC on 
Fjv is constructed with up to 5^-order derivatives with respect to y and t. Using 
the backward-difference method to construct the discretized forms of these derivatives 
requires that we reach back as many as five grid points both temporally and spatially. 
During the initial time-stepping sequence, we only have a history of one time step. 
The current algorithm uses a default of zero if the discretization calls for data at a 
time that is less than zero. This is a poor estimate if the initial condition is some 
non-zero value because it introduces a discontinuity in time. The problem corrects 
itself somewhat if a Sommerfeld condition (J = 1 and Cj = {1}) is utilized on Fjv 
as in Figure 42. However, our goal is to use high-order Higdon conditions to increase 
the accuracy of the estimate. Ebcperimentation showed that a buffer of J grid points 
is necessary to achieve stability for a J^-order Higdon-NRBC. The exact nature of 
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the instability must be furthei investigated. Howevei, for the time being we contdmie 
to present examples while mamtaining the required buffer zone to achieve stability. 
This is one of the reasons Givoh and Neta introduced auxiliary variables [Rjef. 24]. 
See also Givoli et al [Rjef. 26]. 


Truicated DamanD 



Eidendod Danain D 




2 


1 

0 


-1 


u 


-£ 



Figure 42. HNRBC-2DR-2S-1L-U0-V0-T10: Mitigation of Instabhily using Sommer- 
feld Condition (J = 1 and Cj = {!}) 
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C. EXAMPLE TWO: A TWO-DIMENSIONAL MULTI¬ 
LAYER STRATIFIED SCHEME WITH HIGDON NRBC’S 
ON TWO SIDES 


In this section we check the conipatibilil^ of the two-sided Higdon model with 
the multi-layer stratification model presented in Chapter 5. This example will again 
use the semi-infinite quarter-plane (Figure 35) and the same domain specification as 
used in Example 1 of this chapter. The stratified medium is modeled with six layers. 
The layer thicknesses from top to bottom are Oj = {.01, .01, .01, .01, .01, .05}. The 
density for each layer is given hy p. = {1,1.05,1.1, 1.15,1.2,1.25}. A gravitational 
parameter ^ = 10 and a dispersion parameter / = .5 is used. 

The boundary function 17 ^^. on is stipulated to simulate two geophysical 
events. A surface disturbance aldn to the wind is initiated in Li. Using the function: 

h{y,i) = 51 An CCS siii(uj„t), (VI.5) 


(VI.S) 


we let: 

if t < 12.5 Ly< 3.75, 
otherwise 

where Am, Tim, tJm aje given by (VI.3). A second disturbance, simulating seismic 
activity on the ocean floor, is initiated in the Lq using the function: 


lo 


My,i) = .001 cos f-- 1 ^- 


(VI.7) 


Vxv,{y,v = I 


(VI.3) 


and letting: 

if5<t<3^il<^< 2.75, 

0 otherwise 

An other values for 17 ^^., are zero. The simulation is run for 15 time units and results 
are reported using a similar presentation as that used in Example 1. 

At f = 6 (Figure 43), the surface effect has passed through Fjv and F^ and 
||e(f)||n at the surface falls within acceptable tolerances. The effects of the bottom 
disturbance are readily visible near the west boundary. At t = 15 (Figure 44), the 
bottom disturbance has propagated through and has successfully passed through 
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Figure 43. HNREa2DR-2S-6L-U0-V0-T06 with J = 5 and = {1,1,1,1, 1} for 
each Layer 


both artihdal boundaries. Some differences between aJid aje visible, but 

again, the error measure ||e'(t)||j^ is less than the discretization error. From this 
example, we conclude that the multi-layer stratification model will worlc effectively 
when Higdon NELBC’s are placed on two-sides of the domain 
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Figure 44. HNREO2DR-2S-6L-U0-V0-T15 with J = 5 and Cj = {1,1,1,1, 1} for 
each Layer 
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D. EXAMPLE FOUR: A TWO-DIMENSIONAL MULTI¬ 
LAYER STRATIFIED SCHEME WITH HIGDON NRBC’S 
ON FOUR SIDES 


In this example we considei a, truncated domain £7 where Higdon NRBC’s are 
imposed on four sides. The "boundary conditions on Fjv and Fe are discretized as 
before. The boundary conditions on and FiJ^.- are discretized in a similar manner 
(see Section III.C), however forward- vice backward-difference equations are used 
in approximating the spatial derivatives. As before, ^2 is a 5 x 5 square with a 
20 X 20 mesh. Higdon NRBC’s are located at t = 0, ^ = 0, t = 5, and y = S with 
Cj = {1,1,1,1,1} at each boundary. 

The extended domain V is the infinite plane. It is represented by a 15 x 15 
square with a 60 x 60 mesh. The domain of interest £7 is located in the center of 
27 at 5 < J:,y < 10. Higdon boundaries are placed at x = 0, ^ = 0, t = 15, and 
^ = 15 with Cj = {1,1,1,1,1,1,1} at each boundary. This is done for computational 
purposes only and any spurious reflection from these extended boundaries will not 
significantly pollute the domain of interest. On both domains Ax = Ay = .25 and the 
time-step At = .0125. A gravitation parameter of ^ = 10 and a dispersion parameter 
of / = .5 is used. The layer thicknesses are 6^ = (.01, .01, .01, .01, .01, .05} with 
respective densities given hy p. = {1,1.05,1.1,1.15,1.2,1.25}. 

Random values that represent physical disturbances in £7 are introduced at 
designated times at specified x- and ^-ranges in specified layers. At the given time, 
these random values are added to existing values of t/ in the designated layer. The 
dispersive mechanism defined by (II.SS) is than allowed to act on the event. These 
random events are sufficiently complex in nature and “flex^ the Higdon NRRC’s 
exposing any unacceptable spurious reflection. 

For this example two separate disturbances or events are imposed. Elvent 1 is 
a surface event given "by; 


£^.025 


.0001 + rand(—.b, .5) if 1.5 < x, ^ < 3.5, 
0 otherwise 


(VI .9) 
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wheie i 0 pTGsents a disturbance initiated in the iirst layei Lj. at t = .025 

and rand{—.S, .5) is a landoni number on the interval [—.5, .5]. Elvent 2 is a bottom 
event given by: 

r .00015 + rand{-.2S, .75) if 1.5 < t < 2.25 Ll.b<y< 3.5, 

i \ ^ / - - -if - , 

LO otherwise. 

where y) indicates that the event was initiated in at t = 5. The maximum 

amplitude of each event is set to be a small fraction of the total thidoiess of the layer 
in which it is initiated, thus maintaining the integrity of the perturbation analysis 
that was utilized to derive the Klein-Gordon equation. A buffer of at least 5 zero- 
valued grid points was maintained between the NR.BC and the event acknowledging 
the stabdhty considerations discussed in Section VLC, All events must be shifted 5 
units in the x- and ^-directions when initiating activity on V, in order to properly 
place them in the domain’s center. The random seed for these events are preserved 
for use in later examples. 
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Figure 45. EINRBC-2DR-4S-6L-U0-V0-T01: Surface Disturbance Initiated. 
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Figure 46, HNRBC-2DR-4S-6L-UC^V0-T06: Surface Disturbajice has Left Bottom 

Disturbance Initiated. 


A trial is run for 15 time units. At t = 1 (Figure 45), the surface event 
has been initiated and the resulting waves are propagating toward the four sides of 
fl, The error measurement ||e’(t)||n is near zero. At t = 6 (Figure 46), the surface 
event has left £7 with no significant spurious reflection at either the boundaries or the 
corners. Additionally, the bottom event has been initiated and is now propagating 
outward in £7. At t = 13 (Figure 47), both events have passed through the artificial 
boundary successfully. An outward radiating wave front can be seen in the upper- 
right subplot of the euctended domain 17. The left-side sub-plots are remarIcably 
similar indicating that the boundary conditions are suitable for the layered problem. 
At t = 15 (Figure 4S), the wave energy resulting from the disturbances have left 
£7. The effect due to spurious reflection starts to become significant thus increasing 
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Figuie 47. HNRBC-2DR-4S-6L-U0-V0-T13; All Events have Passed thiongh HN- 
EBC’s. 
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Figuie 48. HNRBC'2DR-'4S-6L-U0-VD-T15: Sjmiioiis Reflection Piedominates in Q. 


98 























VII. LINEARIZED SHALLOW WATER 
EQUATIONS WITH NON-ZERO ADVECTION 


In the tMs thaprtei the hneajized sh allow-watei model is fnithei developed hy 
lifting the zeTO-a,dvection leqniienient. RjecaJl the hneajized form of the shallow wa,tei 
model is described by (11.77 and 11.79): 

dkg BHo , Brj' 




= -g 


dz 






dz dz 


Tr + ''ar + ''^ + «''+“) 


+ <™'> 


Bn ^,Bn ,.Bn „ fBu’ Bv'\ 

where U and V aje constant advection terms in the j- arid ^-directions respectively. 
As per discussion in Section II.B.S, if the bottom contour hs is hat, then Hq is 
constant . This allows simpilihcation of (VII.l) to: 


du* 


= -9 


A 

Bz ’ 


Bt 


Bv 


B^ 


, t 


+ + /(f^+ 


3y 


= -9 


&y’ 


(VII.2) 


. v-^ . /5,.* St-’ 


Now define the operator: 


— - Jl— V— 

Dt (9t By' 


= 0 . 


(VII.3) 
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Substntutirig this into ('VII.2) yields; 




= -^ 5 ^’ 




= -g- 


(VIL4) 


D , , „ (dv." dv‘\ 

+ + - 0 , 

wlieie fV and /t/ are constant products of the Coriolis parameter and the advection 
components. Note that this is analogous to the zero-advection form of the shallow 
water equation given hv (II.7S) and (II.SO) with — replacing Equation (VII.4) 
Is combined into a single form using the same steps delineated in Section II.B.9. The 


result is: 




(VII.5) 


where Cq = ^gHQ. This is the Klein-Gordon form of the linearized shallow water 
equation with non-zero advection. This can be rewritten as: 


'^(v) - 1^0 V'l? + fv = S{z,y,t), 


(VIL6) 


sudi that: 


^ 3s „as ,.as „ 


We shall consider (VII.6) in its homogeneous form: 

+ A = 0- 


(VII.7) 


(VII.8) 


Applying the operator listed in (VII.3) twice yields; 

-= — -h - -4- 2U --h 2V--h 2?7V- 

Di^ at^ ^ dx^ ^ dxat ^ dydt ^ dxdy 


(VII.9) 
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Substntutirig this into (VII.S) yields; 

(VII,10) 

This is aji expanded Klein-Gordon equivalent for the linearized shallow watei equa¬ 
tions with non-zeio advection teinis U and V. It applies to a single-1 ayei model, hoit 
will latei be extended to the AT-layer stratified model. 



A, DISCRETIZING THE LINEARIZED SWE WITH CON¬ 
STANT NON-ZERO ADVECTION TERMS 

The following central-diEerence approximations are used to discretize (VII.IO); 






(VILll) 


+ 0(A^^ Ay^). 


Analogous forms are used for approximating and 77 ^, Substituting these 

into (VII.IO) yields; 


+ - + VpJS-l) 


+ - ’/p+u - ’S-u + 


(VIL12) 


+ “ ^P^+1 “ + ^P^-l) 


UV 
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This is the discrete foim of the hneaj shallow wa-tei equation with non-zero advectdon 
terms. Terms that contain an n-j-1 superscript are unknown 17 values which must be 


solved. Making the following substitutions: 


1 

r, 



Ai= ’ 


U 



2 A 1 AV 

2AyAi^ 

F 

2AxAy' 


and moving known terms to the right side of (VII. 12) yields: 


(VII,13) 


+ = 


(2^ + 2B + 2C - - B 

(VII.14) 

“ ^ + L “I" “I" ^ ~ “I" ^ (^^ + 1- “ ■^P-fl-O 

~ ^ ^^+L.^+L ~ ^+L.^-L ~ ^p-L^ + L "I" ^p-L^-l) ■ 

Since there are five unknowns, (VI1,14) must be solved implicitly, A system of N 
equations with N unknowns is set up where AT is the total number of grid points. To 
make the indices compatible with those used in the discretized form of the Higdon 
NRHC (see Sections IV.C and IV.D) we replace n with — 1. This yields: 


(VII.15) 


The Higdon NHBC (III,4) is now applied to the boundaries of the problem. The 
unknowns in the discretized form of Higdon NRBC (IV, 19) are those terms whose 
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□p 0 TatoT Pm does not contain a tdnia shift S^. . For exampla, consider the discretized 
form of a Higdon NRBC with order J = 2 on the eastern boundary (IV.IS): 

[aiaj -|- die's-|- d]_a2S(. diD'^Sf. ^ 

+ die^SfS- + eidsS" + eid^SfS- eie2S-^]rj^^ = 0. 

Shifting the known quantities, e.g. all terms with time shift operators, to the right 
side of (VII. 16) yields: 

[dLds + aie2S~ -\- eia2S~ -\- eie2S~^]7j^^ = 

(VII.17) 

- [dL^isSr + dia2Sr + did2Sf^ -\- die2SrS~ + eid2Sf 
or in alternate form; 

+ (^1^2 + = 

(VII.IS) 

- (dL^^2 -I- dia2)Ti^~^ - (die2 + eid2)Ti^z\^ - did2rf^~J. 

In general, the number of unknowns in a J^-order Higdon NRBC using first-order 
difference approKimations is J -1- 1. The number of unknowns in a J^-order Higdon 
NRHC using second-order difference approximations is 2J -f 1. Equation (VII.15) 
together with an analogous form of (VII.18), adjusted for the Higdon order J, is used 
to set up the system of equations that will be solved to “time-step^ the system. 

B. NUMERICAL EXAMPLE: TWO-DIMENSIONAL SINGLE¬ 
LAYER SCHEME WITH HIGDON NRBC’S ON FOUR 
SIDES WITH NON-ZERO ADVECTION 

In this example, the truncated domain £7 as described in Section VI.D with 
Higdon NRBC’s on four sides is used. As before, the expended domain 27 is an infinite 
plane represented by a 15 X 15 square with a 60 X 60 mesh. £7 is located in the center 
of P at 5 < J,?/ < 10. Higdon boundaries are also imposed on V for computational 
purposes. Spurious reflection from these boundaries should not significantly pollute 


(VII.16) 


103 



£7. On "bat}! domains Ax = Ay = .25 and At = .0125. A gravita-tion pajametei of 
g = 10, dispeision pajametei of / = .5, and a single layer of thickiness 0 = .1 witli 
density y3 = 1 is used. Advection constants of [/ = .5 and V = —.25 are utilized. 

Physical disturbances in are initiated via two separate events. Event 1 is 
given by: 


t=.025 


(VII.19) 


[ .0001 + rand{—.b, .5) if 1.5 < x,^ < 3.5, 

I 0 otherwise 

where S'^'^'^^{x,y) represents a disturbance initiated at t = .025 and rand{—.S, .5) is 

a random number on the interval [—.5, .5]. Event 2 is given by: 

_ r.000015 + rand(-.25,.75) if 1.5 < x < 2.25 1.5 < ^< 3.5, 


lO otherwise. 

^VII.20) 

where S^^^{x,y) represents a disturbance initiated at t = 5 and ra;id(—.25, .75) is 
a random number on the interval [—.25, .75]. Note that these events are analogous 
to the events described in Section VI.D. The random seed was preserved so that the 
same values would be generated for each event. A buffer of at least 5 zero-valued grid 
points was maintained between the NRBC and each event for stability purposes. The 
events are shifted 5 units in the positive x- and ^-directions on 77 in order to properly 
place them in the domain’s center. 

Before running an example, consideration was given to the selection of Cj’s. 
Several experiments were conducted with results reported in Figure 49. Initially, a 
Higdon NREC with order J = 5 and Cj = {Cq, Cq, Cq, Cq, Cq} where Cq = ^5© = 1 
was considered. This was compared to a case where the Cj’s are corrected for advec¬ 
tion. The predominate speed of the gravity wave is Cq. This is affected somewhat 
the dispersion and wave height. However, with the inclusion of advection, the pre¬ 
dominate wave speed with respect to each boundary is aifected more significantly. 
Therefore the C^’s on each boundary are adjusted. These adjustments were made as 
follows: 


= Cj + t/, 

C""*’' = Cj+ V, 


C“' 


Cf 


trfPi 


= Cj - [/, 

= Cj - V. 


(VII.21) 
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Foi this Gxanipl 0 , the adjusted C/s are; 


= {1.5,1.5,1.5,1.5,1.5}, C""*'' = {.75, .75, .75, .75, .75}, 

= {.5, .5, .5, .5, .5}, Cf""* = {1.25,1.25,1.25,1.25,1.25}. 

The results of both runs show asigniiicajit decrease in ||e’(f)|L to about 10“^ atf = 10, 
This eiioi, however, can he reduced further, 




Figure 49, Plot A: J = 5 with Cj = {1,1,1,1,1} adjusted for advection compared 
to J = 5 with Cj = (1,1,1,1,1} unadjusted for advection. Plot B: J = 3 with 
Cj = (1,1, 1} compared to J = 5 with Oj = (1,1,1,1,1}, hoth cases adjusted for 
advection. Plot C; J = 3 with Cj = {.S, .9,1} compared to J = 3 with Cj = {1,1,1}, 
both cases adjusted for advection. Plot D; Corner check for J = 3 with Cj = {.S, .9,1} 
adjusted for advection 


Recall from Section VI,C that a buffer of J zero-vaJued grid points was neces¬ 
sary to achieve stability for a order Higdon NRBC, Wlien advection is incorpo¬ 
rated into the problem, this buffer zone moves horizontally toward at least one of the 
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boiuidaries. TheiefoiG the buffer is compressed with respect to tbe boundary toward 
wbicb it is moving. In order to maintain stability we must eitber increase the size of 
the buffer zone, or reduce the order J. In plot B (top right) of Figure 49 a 5^^-order 
Higdon NEbBC is compared to a 3^^-order Higdon NRBC. In both cases the Cj’s are 
adjusted for advection. In this example, ||e'(t)||j^ is reduced by an order of magnitude 
to about 10““^, 

One further adjustment is possible to reduce ||?(t)||n. Geometric dispersion 
is another factor in the boundaries response to an impinging wave. A wave striking 
normal to the borundary will generally have a wave speed that is approximately Cq. 
In all other cases, the wave speed is less than Cq. An example was set up for J = 3 
in which Cj = {.3, .9,1} with the reduced values taking into account the geometric 
dispersion. Adjusted for advection, the C/s used for the problem are; 


= {1.3,1.4,1.5}, = (.55, .65, .75}, 


= {.3, .4, .5}, = {1.05,1.15,1.25}. 


trfh. 


(VII.22) 


In Plot C (bottom left) of Figure 49 an additional reduction in ||e’(t)||n is evident. 
Further analysis is necessary to determine how to best adjust Cj values for geometric 
dispersion. 

The question of the corner points of Q is again salient in the advection case, 
because the values for Gj on each boundary are now different. Rjecall that there are 
two ways to approximate the boundary values when numerically solving the problem. 
Both approaches are tested here. In the first run the j-boundaries were computed 
first (including the corner points) and the ^-boundaries computed next (excluding 
the corner points). In a second experiment the procedure was reversed and corner 
points were included in the ^-boundaries. Plot D (bottom right) of Figure 49 reveals 
that the solutions are identical. Hence, as concluded earlier, no special handling at 
the corner points is necessary. 

With these results in mind, Higdon NREC’s of order J = 3 with Cj = {.S, .9,1} 
are used. With [/ = .5 and V = —.25, the adjusted C/s are those fisted in (VII.22). 
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A trial is lun foi 10 time units. At i = 1 (Figure 50), event 1 has heen propagating 
outward in for approximately 1 time unit. The effect of advection is apparent 
as the propagation of the gravity wave is tending toward the southeast (i.e. in the 
< .5,—.25 > direction). The leading edge of the wave has passed through the 
but the error measurement is still very small. At t = 2 (Figure 51), event 1 has 
crossed F s aJid Later, at i = 3 (Figure 52), event 1 has crossed Fjv aJid At 
t = 5 (Figure 51), most of event 1 has left FL We note some spurious activity on the 
western boundary. 

At t = 6 (Figure 54), the waves generated by event 2 are approaching and 
Fiv. F^''ent 1 has passed through all four boundaries relatively unperturbed. The plot 
of V reveals that the wave front continues to tend toward the southern and eastern 
portion of the extended domain. 
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Figure 50. HNREG-2DR-4S-1L-U.5-Vm.25-T01; Ev^ent 1 Initiated. 


At i = 10 (Figure 55), the second event has passed through the boundary, The 
wave propagation pattern continues to “drift^ in the direction of advection as revealed 
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Figuis 51, HNREC-2DR-4S-1L-U,5-Vm,25-T02; Ev^ent 1 Crosses Vs aJid T^, 


bv tlie uppei-inght plot of V. Close inspection of the contours reveal spreading where 
the gravity wave is traveling in the direction of advection and compression where the 
gravity wave is traveling against the direction of advection. In the latter case, this 
indicates a steeper wave front. Since the gravity wave is omni-directional, this effect 
varies throughout the plot, In the noise of spurious reflection is now visdhle. 

This experiment was repeated for two other sets of values for U and V. In 
the hr St variation (Figure 56) the magnitude of the advection constants were lowered 
to U = A and V = —.15. As expected, there is a decreased tendency toward the 
southeast. Also notable is a reduction in the error measurement. In the second 
variation (Figure 57), the magnitude of the advection constants was increased to [/ = 
.6 and V = —.35, The tendency to the southeast, as well as the error measurement, 
has increased. These results indicate that the model is behaving as expected with 
regards to the rate and direction of advection. However, as the magnitude of the 
advection constants is increased, the measured error will also increase. In the current 
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Figuis 52. HNREC>2Da-4S-lL-U.5-Vm,25-T03; Ev^ent 1 CiQsses Tjv and 


exampile, the magnitude of the advection is 4 to 7 times greater than the magnitude 
□f the depth. In a real world problem, where the open ocean is the medium of 
propagation, advection constants are expected to be significantly smaller. 
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Figuie 53, HNRJBC-2DRr4S-lL~U.5-Vm.25-T05; Rv&nt 1 Lea^’es Q with Visible Spu- 
lious Rejection at F, 
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Figure 54, HNaBC'2DR-4S-lL-U.5-Vm,25.T06; Event 2 Initiated. 
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FiguTG 55 . HNREa 2 DR- 4 S-lL-U. 5 -Vm. 25 -T 10 ; The Noise of Spuiious Reflectdon 
Evident at the Bottom Left Plot. 
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Fi^ie 56 . HNRBO 2 DR- 4 S-lL-U. 4 -Vm.l 5 -T 10 ; End of Rim 
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Figui& 57. HNRBa2D-4S-lL-U.6-Vm.35-T10; End of Run 


112 



















C. LINEARIZED SWE WITH CONSTANT NON-ZERO 
ADVECTION TERMS EXTENDED TO THE N-LAYER 
MODEL 

Recall the linearized form of the shallow water equation with non-zero advec- 
tion terms for a AT-layer stratified model (V.IO): 


I TJ I I.- ffir \ \ — rr f Y" I V" 1 




— +1/ — + v'^ + f(u- +ii') = -o + v'^] 


/s«r stn _ 

^ + t/,^+V,—+ 0,1 —+ —1 -0. 


^th 


where and V] are the x- and ^-components of advection in the -layer. The same 
techniques as those described in Section II.B,8 are used to generate a single equation 
for 7]ij the perturbation of the 7^-layer. As an intermediate step a result that is 
analogous to (V.19) and (VII.6) is generated; 


D' 


Dt^ 




{7^i)-gQiV^ \ ^ 17 ; + ) + /V = Si{x,y,t), 

y=l y=i j 


JV 


where S^ix^y^i) is the source function for the z"^-layer of an iV-layer model. The 
is defined as; 


(VII.23) 


operator 


Dt. 


Siix^y^t) must satisfy; 

D 


Dt. 


^ TT ^ T- ^ 


{Si{x,y,i)) - 0 ^ +^i +'^i -^■ 


(VIL24) 


(YI.25) 
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Consdd 0 Tiiig th 0 hoiiiog 0 ii 0 oiis foim foi (VII.23) ajid 0 xpanding yi 0 lds: 



This is th 0 liiL 0 ariz 0 d shallow wat 0 T 0 q'ua,tioii foi aji TV-lay 01 iiiQd 0 l with 
adv 0 ctioii t 0 iins. 


(VII.26) 


iiQn-z0io 
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D. DISCRETIZING THE LINEARIZED SWE A-LAYER 
STRATIFIED MODEL WITH CONSTANT NON-ZERO 
ADVECTION TERMS 

Astefoie, using the ceutral-diffeience appioodmataous in (VII.11) to discietize 
(VII.26) yields; 




+ 


9 AjAt ^.P-L^ + ^^.P-L^.J 


_L _ Yi _ ) 

^ Vi.pj}+1 Vi.p.i}-1 ^ Vi.pjj-lJ 


UM . 




itr 


E ^ (’7"p+i., - 

j=l ri 


(VIL27) 






^-L 


\j=ip^ 




= 0 . 


This is the discretized form of the 7V-layei hneazized shallow wa,tei equation with 
constant, non-zero advection for the layer. All terms that contain the n H- 1 
superscript are unknown values of 77 ^ which must he solved. Making the following 
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sulDstitutions; 


D- = ^— 

' 2 AxAt 


A = 


1 

At=’ 



Ei = 


v; 

2 Aj(At’ 


Fi 


2ATAt»’ ' Ur;) ’ ’ iAt,J ’ 


(VII.2g) 


maving known teinis to the light side of (VII.27), and shifting the time index yields: 




( 2 ^ + 2B, + 2 C; - 


+ A (C+M - C-iO + A fc'+i - ’?"p1-i) 


- A {nU 


P+L^+L 


V^.pfL^-L 


^.p—L^+L ~l” ^^.p—L^ — 




^-L 


j=L 


+ G', 


i: (’)f.pt^ - + ^7^;^ 


p-l-fl) 


\J=t 


A Hi 


£f teL-2^;.«uci-0 

j = L 


+ A 


L (’7Zp^+i - 


J=^ 


9rf"-^ 
^ h-T^ 



(VII.29) 


As in the single-layei advection case, (VII,29) must he solved imphdtly foi each layei 
Li. The system of equations is completed on the boundaries using the discretized 
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Higdon boimdary equations as discussed in section VII.B. 


E. NUMERICAL EXAMPLE: A TWO-LAYER SCHEME 
INCORPORATING ADVECTION 

In this exajnple the domains and V as described in Section VII .B ■with identi- 
caJ positioning of the Higdon NRBC’s are "utilized, The following pioblem pazameters 
aze used; 

Ax = Ay = .25, At = .1, 

5=10, / = .5, 

0i = {.O3,O7}, Pi = {1,1.05}, 

Ui = {.025, .025}, Vi = {-.025, -.025}, 

J=5, Cj = {.S,.7,.3,.9,l}. 

Collecting "the GjS for advection yields: 

= {.325, .725, .S25, .925,1.025}, = {.575, .675, .775, .875, .975}, 

= {.575, .375, .775, .375, .975}, = {.625, .725, .825, .925,1.025}. 

A single physical disturbance is initiated in fl and is given 1^: 

f .000001 + randi—.S, .5) if 2 < t, v < 3, 

S^"'" = J _ (VIL30) 

10 otheiwise , 

wlieie lepiesents a disturbance initiated in at t = .1 and Tand{—.b, .5) 

is a landoni nunibei on tbe interval [—.5, .5]. The example is run for £ve time steps. 

At t = 1 (Figure 5S), the disturbance has been underway for approximately 
one second. Minimal spurious reflection occurs at the boundaries. In the lower-right 
plot two additional measurements are noted. The first, “Max Ref SurP is 
measured over the entire run. The nex±, “Max He'll Ratio^ is given by: 

Max ||e||Ratio = (VIL31) 

Since both are maximums ex±racted from the data generated over the entire run, 
they will not change with time. Both are used in the next section to compare Higdon 
boundary effectiveness. 


117 



At t = 5 (Figure 59), most of the wave action has left £7. The lesidnal action in 
the truncated domains are, for the most part, similar. There is, however, some visible 
difference near the boundaries resultant from spurious reflection. The lower-right plot 
reports: 

Mas ||e||Ratio = 1.08%. 

That is to say, the maximum error norm ||e|| at t = 5 was 1.08% of the |i7|ma^. As we 
shall see in the next section, this is a favorable measurement for the Higdon NRBG. 
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Figure 58. HNRBc:>2DR-4S-2L-Up025-Vmp025-T01; Disturbance Initiated. 
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Figure 59. HNREa2DR-4S-2L-Up025-Vmp025-T05; End of Run, Some Noise at 
Boundaries of Bottom Left Plot. 
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F. VARYING PARAMETERS FOR THE TWO-LAYER 
SCHEME WITH ADVECTION 

Thus fax in this dissartatioii, the truncated domain has been numeric ally 
compared to its infinite analog, the extended domain D. Unfortunately, most of the 
computational assets are invested in generating a solution for this “pseudo-infinite^ 
domain. In the four-sided Higdon NHBC’s, the extended domain V contained nine 
times as many grid points as Addition of advection considerations required em¬ 
ployment of implicit methods. In these cases, computer storage became alcey factor. 
For example, V for a four-sided problem with AxAy = .25 required a 61^ x 61^ ma¬ 
trix to generate a solution whereas required a much smaller matrix of dimensions 
21^ X 21^. Obviously sparse matrix procedures would alleviate this requirement, but 
the point illustrates that maintaining an extended domain becomes unwieldily if the 
domain of interest is enlarged or the grid is refined. 

In the following subsections, parameters of a two-layer advection problem are 
varied and the resulting measurement of “Max ||e| | Ratio^ (VII.31) is compared. From 
these comparisons, conclusions are listed that predict parameter limits that destabilize 
the Higdon NRHC. Knowledge of this behavior will allow us to drop extended domain 
comparisons and focus attention on the domain of interest. 

In ail triais, the domain is a 5 X 5 square. The grid is varied as specified in 

each set of trials. The foil owing event is used to initiate action in ail trials; 

f .000001 + .5, .5) if 2 < t, ^ < 3, 

^ (VIL32) 

L 0 otherwise , 

where is an event generated in the i^-layer at t = At (i.e. the first time 

increment). Each trial is run for 5 time units and all C^’s are adjusted for advection. 
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1, Varying C'/s, Ajt, and Ay 

In this set of trials the Higdon coefficients Cj were varied. The following 
parameters were fixed: 

= ^, = { 1 , 1 . 05 }, 

Ui = {.025, .025}, v; = (-.025, -.025}, (VII,33) 

At = .1 Ai = Ay = .5 

The value for |i 7 |fnaa; was determined to be 1.92 X 10“^ and was imaffected by the 
choice of C/s. Rjesults for Az = Ay = .5 are posted in Table IV and were compared 
using ^ measure described in Section VILE. This set of numerical trials was 


Table IV, Varying C^-s with Az = Ay = .5 





I^TI mo.! 


3.47% 

{0,,5,1} 

3.83% 


3.53% 

{.7,,8,.9} 

3.97% 

1—1 

no 

3.57% 

{,8, .9,1} 

4.13% 

{,1,75,1} 

3.58% 

{111} 

4.21% 

{.025, .6,975} 

3.58% 

{1,1} 

4.23% 

1—1 

3.61% 

{,1,.1,1} 

5.99% 

{.1,.5,.9} 

3.63% 

{1,1,1,!} 

92.2% 

{,5, .5,,5} 

3.76% 




repeated for Ax = Ay = .25. In this case was determined to be 1.41 x 10“*^, 

a change that reflects the nature of the random event (VII.32) that was used to 
generate the disturbance. The results of these trials are reported in Table V. Several 
conclusions are drawn from these; 

• Refining the grid results in a more efficient boundary conditions (e.g. less 
spurious reflections), 

• Higher order Higdon NRBC’s tended to generate more effective boundary con¬ 
ditions. 

• A buffer zone as discussed in Section VI.B and VII.B must be considered. 
This limits the Higdon order J that can be used. For Ax = Ay = .5, the trial 
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T&,ble V, Vajying C^-s witli Ai = Ay = .25 



■IBI^ 



{,6, .7,, 3, .9,1} 

1.03% 

{.1,.75,1} 

2.14% 

{.9,.925,.95,.975,1} 

1.13% 

{.5,.5,.5} 

2.15% 

{1,1,1,1,1} 

1.16% 

{.1,.6,1} 

2.24% 

{.1,.5,.S„7,1} 

1.19% 

{.1,.5,1} 

2.37% 

{.7,.3,.9,1} 

1.21% 

{.1,.5,.9} 

2.37% 

{1,1,1,!} 

1.24% 

{1,1} 

2.37% 

{.3,. 9,1} 

1.57% 

{.025,.6,.975} 

2.43% 

1 - 1 

1 - 1 

1 - 1 

1.53% 

{0,.5,1} 

2.64% 

{.7,.3,.9} 

1.60% 

{.1,.1,1} 

3.71% 

{.6,. 3,1} 

1.63% 

1 - 1 

1 - 1 

1 - 1 

1 - 1 

1 - 1 

1 - 1 

393% 

{.4,. 7,1} 

1.77% 




became imstable for J = 4. For At = = .25, instability Qccnied at J = 6. 

In the latter case, higher order Higdon NRBC’s could be used increasing the 
effectiveness of the boundary condition. AnxOiary variables (see [Ref. 26]) will 
alleviate this problem. 

• Distributing the Cj’s on the interval [.5,1] seemed to reduce spurious reflection, 
but how to distribute these values best could not be determined. In any event, 
it appeared that at least one of the values should be equal to the gravity wave 
speed Co, which in this case was 1. 

• As approached 5%, spurious reflections became visible. Values of 3% 

produced acceptable results. Obviously, less is better. 

2, Varying U and V 

In this set of trials U and V are varied, however, their respective values are 
kept the same in each layer. The following parameters are fixed: 

J/, = {.03,.07}, ^ = {1.,1.05}, 

At = .1, At = A^ = .5, (VII.34) 

J=3, Cj = {.3,.9,1}. 

Results are posted in Table VI. These trials were repeated for At = A^ = .25 with 
results reported in Table VII. Another set of trials were conducted for Cj = {.1, .6,1} 
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Table VI. Varying V and V with Ax = = .5 ([/, V Equal in each Layer) 


V 

V 

17 

1 ^1 mo. 3- 

{. 0001 ,, 0001 } 

{-, 0001 ,-. 0001 } 

1.92 X 10 -'= 

4.13% 

{, 01 ,, 01 } 

{-. 01 ,-. 01 } 

1.92 X 10-'= 

4,13% 

{,025,,025} 

{-.025,-.025} 

1.92 X 10-'= 

4,13% 

{,03,,03} 

{-.03,-.03} 

1.92 X 10-'= 

4,14% 

{,035,,035} 

{-.035,-.035} 

1.92 X 10-“ 

4,17% 

{,04,, 04} 

{-.04,-.04} 

1.92 X 10-^ 

4,20% 

{,05,,05} 

{-.05,-.05} 

1.92 X 10-'= 

4,31% 

{. 1 „ 1 } 

{-, 1 ,-.!} 

1.93 X 10-^ 

7,01% 

{,25,,25} 

{-.25,-.25} 

2.32 X 10-= 

95.5% 


Table VII. Varying V and V with Ax = iS.y = .25 ([/, V Equal in each Layer) 


V 

V 

17 

1 1^1 1 rTU.3 

Ifjlmai 

{. 0001 ,, 0001 } 

{-, 0001 ,-. 0001 } 

1.42 X 10-“ 

2.28% 

{, 01 ,, 01 } 

{-. 01 ,-. 01 } 

1.42 X 10-^ 

2,24% 

{. 1 „ 1 } 

{-,l,-.l} 

1.40 X 10-^ 

3,98% 

{.25,.25} 

{-.25,-.25} 

1.07 X 10-= 

58.5% 


with Ax = ISy = .25. All other parameters were unchanged. Again U and V was 
varied, but this time the advectdon coefficients were allowed to differ between the 
layers. Table VIII reports the results when the disturbance is initiated in while 
Table IX reports results when the disturbance is initiated in Lj. From these results, 
it can be concluded; 

• The problem became unstable when [I, V became large. A good rule of thumb 
is to keep [/, V (J, the total depth of the medium. This is a physical con¬ 
straint in an ocean environment that assumes that the magnitude of advectdon 
will be much less than the magnitude of depth. 

• Refining the grid improved performance. 

• The instabihly of the problem intensifies when t/, V differs between layers. 
Such problems should be avoided, when advectdve differences between layers 
are small. 
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Tabde VIII. Vajying U and V with Ax = Ay = .25 for Event Initiated in Li {U, V 
not Necessarily Equal in each Layer) 


u 

V 


k\n.^ 

{.02,.02} 

{-.02,-.02} 

1.92 X 10^ 

3.50% 

{0,0} 

Lo2,-.02} 

1.92 X 10-“ 

3.51% 

{.02,.02} 

{0,0} 

1.92 X 10-® 

3.60% 

{0,0} 

{0,0} 

1.92 X 10^ 

3.92% 

{.01,0} 

{-,01,0} 

1.92 X 10-^ 

7.26% 

{0,.02} 

{0,-,02} 

1.92 X 10"“ 

12,6% 

{.01,-,01} 

{-.01,.01} 

1.94 X 10-® 

12,8% 

{.02,0} 

{-,02,0} 

1.92 X 10^ 

13,1% 

{.02,-,02} 

{-.02,.02} 

1.99 X 10-® 

25,5% 


3. Varying Layer Thicknesses in a 2-Layer Problem 

In this set of trials, the layer thicknesses are varied. The following parameters 


were fixed: 


At = .1 


J = 3 

U = {.025,.025}, 


Ax = Ay = .5 
Cj = {.8,.9,l} 

V= {-.025,-.025}. 


(VII.35) 


Tatle DC, Varying U and V with Ai = Ay = .25 foi Event Initiated in L 2 (V, V not 
Necessarily Equal in each Layer) 


u 

V 



{0,0} 

{-,02,-.02} 

1.92 X 10-“ 

3.61% 

{.01,-.01} 

{-.01,.01} 

1.94 X 10^ 

3.63% 

{0,.02} 

{0,-,02} 

1.92 X 10-^ 

3.71% 

{.02,.02} 

{-,02,-.02} 

1.92 X 10-® 

3.74% 

{0,0} 

{0,0} 

1.92 X 10^ 

4.38% 

{.02,-,02} 

{-.02,.02} 

1.99 X 10^ 

6.31% 

{.02,.02} 

{0,0} 

1.92 X 10"® 

6.32% 

{.01,0} 

{-.01,0} 

1.92 X 10-® 

6.85% 

{.02,0} 

{-,02,0} 

1.92 X 10-^ 

11,3% 
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Results foi a distuibance initiated in Li and Lj are posted in Table X. It is readily 


Table X. Varying Layer Thidoiess in a 2-Layer Problem with an Li Initiated Event 
(left) and L 2 Initiated Elvent (right) 


fli 

^ max 

1 kll rminr 



Ifjlmai 


{.1,0} 

1.39 X 10 -^ 

3.57% 

1.39 X 10-'=' 

3.57% 


1.99 X 10-® 

3.67% 

1.39 X 10-'^ 

3.57% 

{.2,.3} 

1.99 X 10-® 

3.51% 

1.39 X 10-'= 

3.53% 

{.3,.7} 

1.92 X 10-“ 

3.47% 

1.39 X 10-“ 

3.53% 

{.4,.6} 

1.92 X 10-“ 

3.45% 

1.89 X 10-“ 

3.53% 

{.5,.5} 

1.91 X 10-® 

3.46% 

1.39 X 10-'= 

3.53% 

{.6,.4} 

1.91 X 10-^ 

3.43% 

1.39 X 10-“ 

3.53% 

{.7,.3} 

1.90 X 10-^ 

3.51% 

1.89 X 10-“ 

3.53% 

{.3,.2} 

1.90 X 10-“ 

3.53% 

1.39 X 10-“ 

3.57% 


1.39 X 10“® 

3.56% 

1.89 X 10-“ 

3.57% 


concluded that; 


• Varying layer thicknesses does not alter the effectiveness of the Higdon NRBC. 
This conclusion is independent of the layer in which the disturbance was ini¬ 
tiated. 


4. Varying Density Distribution in a 2-Layer Problem 

In this set of trials, the density distributions are varied. The following param¬ 
eters were fixed; 

= {.03, .07}, 

At = .1 At = iS.y = .5 

(VII.36) 

J = 3 Cj = {.S,.9,l} 

V = (.025, .025}, V = {-.025, -.025}. 


Results for a disturbance initiated in Z-l and Z -2 are posted in Table XL From these 
trials it is concluded that; 


• Density changes do, to some extent alter the behavior of the Higdon HNRBC, 
but in general, they remain effective as long as the densily increases mono ton¬ 
ic ally with each increasing layer. 
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Tabde XI. Vajying Layer Thickness in a 2-Layer Problem with a Initiated Event 
(left) and Lj Initiated Event (right) 


1 ^ Imaa: 

IfjIniLiT 

1.39 X 10-^ 

3.53% 

1.39 X 10-“ 

3.53% 

1.33 X 10-® 

3.60% 

1.33 X 10-® 

3.60% 

1.33 X 10-^ 

3.57% 

1.91 X 10"“ 

3.47% 

2.10 X 10-“ 

3.20% 

1.39 X 10“^ 

3.92% 

1.93 X 10-^ 

3.92% 

1.09 X 10-^ 

6.94% 


pi 

^ max 

11^ 11 ma^ 

IfjImoT 

{1,1} 

1.39 X 10-^ 

3.53% 

{1,1.05} 

1.92 X 10-“ 

3.47% 

{1,1.10} 

1.99 X 10-'= 

3.20% 

{1,1.15} 

2.03 X 10-'= 

2.02% 

{1,1.20} 

2.05 X 10-'= 

2.93% 

{1,1.25} 

2.05 X 10-“ 

3.02% 

{1,1.50} 

2.13 X 10-“ 

3.11% 

{1.05,1} 

1.35 X 10“^ 

4.35% 

{1.10,1} 

1.92 X 10-^ 

3.34% 

{1.25,1} 

2.10 X 10-* 

3.69% 


5, Varying At in a 2-Layer Problem 

In this set of trials At is varied. The following parameters were fixed; 


= {.03,.07}, ^ = {1,1.05} 
Ax = Ay = .5, 

J = 2, C^ = {.3,.9,1}, 

y = {0,0}, v- = {0,0}. 


(VII.37) 


Results are posted in Table XII. In general, decreasing At does not change the 


Table XII. Varying At in a 2-Layer Problem with a L]_ Initiated Event 


At 


11 ^11 mas‘ 

1^1 mar 

.1 

1.92 X 10-“ 

4.13% 

,05 

3.32 X 10^ 

3.13% 

.01 

1.91 X 10-= 

3.25% 

.005 

3.32 X 10-= 

3.23% 

.002 

9.57 X 10-= 

3.23% 


effectiveness of the Higdon NRBC. However it appears that liylmoE is inversely pro¬ 
portional to At. The probable reason for this behavior is that the introduction of 
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the event pi□ duces a discontdnnity in tj with regards to A smaller At results in a 
more pronounced discontannitv. It takes several numerical time-steps to sort out the 
discontinuity. The result is a that is greater then the maximum amplitude of 

the event induced on Q. 

G, THE “HIGDON MATRIX” 


hWm Irr^-i farJl «li 



s iDD I9D aiD an no an m 


Figure 60, Higdon Matrix Image for (20 x 20) with Higdon NRBC’s with Order 
J = 9 Applied to Four Sides 

As mentioned in Section VILA, when non-zero advection terms are incorpo¬ 
rated into the shallow water model the problem must be solved implicitly, and a 
X matrix with a bandwidth of 2N^ is generated. An image of this matrix 
is presented in Figure 60 where zero elements are black and non-zero elements are 
white, Here the truncated domain ^2 is approximated using 21 x 21 grid and Higdon 
NRBC’s of order J = 9 are applied to all four sides. On the top and the bottom of 
the image we see 10 light diagonal lines. These lines represent the discretization for 
the ^-boundaries Fjv (top) and Fjv (bottom). The heavier line along the diagonal is 
three points thick and is flanked to the left and right by two thinner lines. These 
result from the discretization of the interior points. Finally, the periodic “ishort-spikie^ 
pointing to the left and right were generated by the Higdon NRHC’s on F^ and FiJ^/ 
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i0spGctdvGly. NotG that thGIG aiG only 19 Gach of thGSG shoit horizontal hnGS. This 
indicatGS that thG coinGi points wgig inclndGd in thG ^-bonndariGS, othGiwisG 21 (N^) 
of such pairs would bG visiblG. 

It is GvidGnt from thG imagG, that thG Higdon matrix: rGquiiGd for non-ZGio 
advGction problGmis sparsG. ThG numbGi ofnon-ZGio points gGUGiatGd by thG domain 
intGrior is; 

5(/\r,-2)(«,-2), (VIL3S) 

and thG numbGi of non-ZGio points gGUGiatGd by thG four Higdon NHBC’s is; 

(2N^ + 2Ny - 4)( J + 1). (VII.39) 


ThGiGfoiG thG fraction of non-ZGio Glcmcnts in the matrix is; 


5 (Ar, _ 2)(Ny - 2) + {2N^ + 2Ny - 4)( J + 1) 


(VII.40) 


In thG casG of our GxamplG whGiG = Ny = 21 and J = 9, only 1.34% of thG matrix 
is populatGd with non-ZGio yaluGS. In thG AT-layGi modGl, N of such matricGS might 
bG produced severely taxing computer memory. Increasing domain size and a refining 
grid would further exacerbate the problem, Clearly sparse matrix: procedures are in 
order and should be investigated to stream hne the non-zero advection scheme and 
produce a faster algorithm. 
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VIII. CONCLUSION AND 
RECOMMENDATIONS FOR FURTHER 

RESEARCH 

Fiom the pieceeding investigataon we conclude that a very large domain with 
dispersive wave action governed by linear SWE’s can be effectively restricted using 
Higdon NRBC’s. This boundary condition will also worlc for versions of the linear 
SWE that include the effects of stratification and advection. However, as the SW’E 
model becomes more complex, careful consideration must be given to values assigned 
to problem parameters to ensure the stability of the scheme. 

There are many aspects of the problem addressed in this dissertation that 
should be investigated further. With regards to the SWE and the geophysical envi¬ 
ronment, further model development should include: 

• Discretization of the non-linear SWE. 

• Incorporation of a non-constant Coriolis parameter. 

• Inclusion of a bottom topography that varies with x and y. 

• Consideration of terms resulting from the Earth’s curvature (i.e. domains with 
horizontal dimensions greater than 1000 km). 

With regards to computational techniques, one should consider using sparse matrix 
algorithms for models requiring implicit solution techniques. Finally, with respect to 
the Higdon NRJBC, the following areas of research are recommended: 

• Development of schemes that use auxiliary variables to eliminate non-zero 
buffer zones, This is done outside this dissertation, 

• Development of schemes to optimize selection of Higdon coefficients. 

• D e vel opment of meth odstoupdateHigdonco efficients dyn ami c ally / ad ap tively. 

• Extension of the rectangular domain to three dimensions. 

• Development of Higdon NRBCs for cylindrical and spherical coordinates. 
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• Cieataon of two-way Higdon NRJBC’s to incoipoiate data fiom the immediate 
vicinily of the trimcated domain. 

These extensions to the cniient research are potentially of great value to oceanogra¬ 
phers and meteorologists alike and may further promote the use of Higdon NEtBC’s 
in weather prediction models. 
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